Минимальный базис задач для валидации методов численного моделирования турбулентных течений вязкой несжимаемой жидкости Текст научной статьи по специальности «Физика»

МЕХАНИКА ЖИДКОСТИ, ГАЗА И ПЛАЗМЫ

УДК 532.54

А.С Козелков1'2, Ю.Н. Дерюгин1, Ю.А. Циберева1, А.В. Корнев3, О.В. Денисова1, Д.Ю. Стрелец3, А.А. Куркин2, В.В. Курулин1, И.Л. Шарипова1, Д.П. Рубцова1, М.А. Легчанов2, Е.С. Тятюшкина1, С.В. Лашкин1, А.В. Ялозо1, С.В. Яцевич1, Н.В. Тарасова1, Р.Р. Гинниятуллин1, М.А. Сизова1, О.Л. Крутякова1

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

ФГУП «РФЯЦ-ВНИИЭФ» 1, Нижегородский государственный технический университет им. Р.Е. Алексеева ,

ОАО «Компания Сухой»

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

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

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

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

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

Введение

Быстрое развитие в последние десятилетия различных областей науки и техники отчасти явилось следствием колоссального прогресса в применении численных методов в исследовании гидродинамики жидкостей, как теоретической, так и экспериментальной. Вычислительная гидродинамика по своей полноте, разнообразию задач, и что особенно актуально, практическому применению выделилась в отдельную самостоятельную науку. Практическому внедрению вычислительной гидродинамики (CFD-Computational Fluid Dynamics) способствует бурный рост производительности вычислительной техники, введение в строй суперкомпьютеров терафлопного и даже петафлопного класса, а также серьезные успехи, достигнутые в последние годы в области построения эффективных численных методов, которые позволяют успешно решать достаточно широкий класс практически важных задач в различных прикладных областях.

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

© Козелков А.С, Дерюгин Ю.Н., Циберева Ю.А., Корнев А.В., Денисова О.В., Стрелец Д.Ю., Куркин А.А., Курулин В.В., Шарипова И.Л., Рубцова Д.П., Легчанов М.А., Тятюшкина E.C., Лашкин С.В., Ялозо А.В., Яцевич С.В., Тарасова Н.В., Гинниятуллин Р.Р., Сизова М.А., Крутякова О.Л., 2014.

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

• необходимостью снижения сроков разработки и стоимости новых образцов;

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

Рис. 1. Основные задачи при моделировании систем самолета

Рис. 2. Классы практических задач и соответствующие задачи валидации

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

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

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

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

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

Изложение основных принципов верификации и валидации методов CFD-моделирования представлены в [1, 37].

Процедура валидации для CFD пакетов включает в себя несколько этапов:

1. Оценка сходимости итераций. При численном моделировании выделяют два типа сходимости - сходимость по итерациям (по времени) и сходимость по пространству. Большинство практических задач решаются в «квазистационарном» режиме, когда моделирование осуществляется до того момента, когда решение устанавливается и больше не меняется со временем (на последующих итерациях). Такую сходимость, называемую сходимостью по итерациям (времени), определяют с помощью достижения определённого значения невязок численного решения. Если учитывается реальное физическое время процесса, то сходимость по времени определяется путем оценки достижения численным решением конечного времени на каждом временном шаге. Сходимость по пространству определяет погрешность дискретизации исходной системы уравнений в физико-математической модели. Оценка сходимости по пространству осуществляется выполнением моделирования на нескольких последовательно сгущающихся сетках. По мере измельчения расчетной сетки погрешности дискретизации по времени и пространству асимптотически должны приближаться к нулю.

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

3. Сравнение с экспериментом. Самый важный этап валидации. Экспериментальные

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

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

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

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

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

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

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

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

1. Оценка погрешности расчета

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

- р I

§ = \ ан / эксп_РаА . 100%, (1)

ан/эксп

где Ган/эксп - аналитическое или экспериментальное значение; ¥расч - расчетное значение исследуемой величины.

Для оценки абсолютного отклонения измеренных значений от среднеарифметического используют величину среднеквадратического отклонения, вычисляемого по формуле[2]:

У (р - р

/ у V ан / эксп ра

\ 2

Л (N -1)N (2) здесь N - количество проведенных численных экспериментов.

2. Базис задач валидации для течений вязкой несжимаемой жидкости

Течение вязкой несжимаемой теплопроводной жидкости описывается физико-математической моделью, основанной на системе уравнений Навье-Стокса [3] и характеризуется следующими физическими свойствами:

1) течение может быть как стационарно, так и нестационарно;

2) течение может быть как ламинарным, так и турбулентным;

3) течение содержит отрывные зоны и вторичные токи;

4) течение переносит тепло и участвует в теплообмене с твердым телом;

5) в течении может преобладать как вынужденная, так и естественная конвекция.

Минимальный валидационный базис должен содержать задачи, на основании решения

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

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

Таблица 1

Таблица задач валидационного базиса

№ Название задачи Учитываемые свойства Оцениваемые параметры Доступные данные

1 2 3 4 5

Ламинарные течения

1 Ламинарное течение Стационарность, лами- Профиль Аналитическое

между пластинами нарный режим течения скорости решение

2 Ламинарное течение Стационарность, лами- Перепад давления, Аналитическое

в трубе треугольного сечения нарный режим течения профиль скорости решение

3 Движение жидкости в квадратной каверне с движущейся стенкой Стационарность, лами- Профиль Аналитическое

нарный режим течения скорости решение

Окончание табл. 1

1 2 3 4 5

4 Безотрывное обтекание шара, задача Стокса Стационарность, ламинарный, режим течения Сила сопротивления Аналитическое решение

Турбулентные течения

5 Течение в круглой трубе Стационарность, переходный и турбулентный режим течения Перепад давления, профиль скорости Аналитическое решение

6 Турбулентное обтекание пластины Стационарность, турбулентность Коэффициент трения Эксперимент

7 Течение в круглом канале с внезапным расширением Стационарность, турбулентность Перепад давления Эксперимент

8 Течение в плоском асимметричном диффузоре Стационарность, турбулентность, отрыв Профили скорости Эксперимент

9 Течение в отводах круглого и квадратного сечений Стационарность, турбулентность Перепад давления Аналитическое решение

Нестационарные течения

10 Развитие ламинарного течения в трубе Нестационарное течение, ламинарный режим течения Профили скорости Аналитическое решение

11 Ламинарное обтекание цилиндра Нестационарное течение, ламинарный режим течения Число Струхаля, лобовое сопротивление Эксперимент

12 Поперечное обтекание трубного пучка Нестационарное течение, турбулентность Профили скорости Эксперимент

13 Вырождение однородной изотропной турбулентности Нестационарное течение, турбулентность Энергетические спектры Эксперимент

14 Развитое турбулентное течение в канале Нестационарное течение, турбулентность Профиль скорости, напряжения Рей-нольдса Прямое численное моделирование

Смешанная и естественная конвекция

15 Конвекция в плоском канале Смешанная конвекция Профили скорости Аналитическое решение

16 Конвекция в квадратной полости Естественная конвекция Профили скорости, число Нуссельта Эксперимент

17 Развитая турбулентная конвекция Турбулентная естественная конвекция Профили скорости Эксперимент

Течения с теплообменом

18 Течение в канале с обратным уступом с подогреваемой стенкой Стационарность, турбулентность, отрыв, теплообмен Коэффициент трения, число Стэнтона Эксперимент

19 Охлаждение твердотельных блоков в плоском канале Сопряженный теплообмен, объемное энерговыделение Профиль температуры Кросс-верификация

20 Вентилируемая квадратная каверна с твердотельным блоком Сопряженный теплообмен Изотермы температуры, число Нуссельта Кросс-верификация

3. Описание задач базиса валидации

Приводится детальное описание задач базиса валидации, представленных ранее. Также здесь приводятся результаты валидации пакета программ ЛОГОС для данного класса задач. Пакет программ ЛОГОС — отечественный пакет программ инженерного анализа, предназначенный для решения сопряженных трехмерных задач конвективного тепломассопере-носа, аэродинамики, гидродинамики и прочности на высокопараллельных ЭВМ [4-5].

Ламинарные течения

Задача 1. Ламинарное течение между пластинами

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

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

Таблица 2

Числа Рейнольдса

АР -3 -2 -1 0 1 2 3

Re 1908 1240 570 100 770 1440 2108

Расчетная сетка содержит 250*20*1=5000 ячеек. На входе и выходе задается перепад давления - от Р1 = -3Па до Р2 = 3Па (рис. 3). На нижней и верхней стенке задано ГУ «Стенка без проскальзывания», на верхней стенке также задана скорость. На передней и задней поверхностях расчетной области задается условие симметрии.

Рис. 3. Геометрия задачи

Для оценки правильности численного решения необходимо найти профиль скорости в канале, для которого существует точное аналитическое решение [6, 7]:

и _у И ■ У Ф(1_у) (3)

и И 2и■у йх И Кривые распределения скоростей, даваемые решением уравнения (3) для различных значений перепада давления, определяют тип течения.

Форма кривой распределения скоростей определяется безразмерным градиентом давления [6, 7]:

Р = _И!_ (_ ± 1. (4)

2ц,и ^ йх у

Для Р > 0, т.е. для случая перепада давления в направлении движения верхней стенки, скорость положительна по всей ширине канала. При отрицательных Р в некоторой части поперечного сечения возможны отрицательные скорости - возвратное течение. Вблизи неподвижной стенки такое течение возникает уже при Р < — 1. Объясняется это тем, что для частиц жидкости, находящихся вблизи неподвижной стенки, увлекающее действие соседних более быстрых слоев не в состоянии преодолеть перепад давления, действующий в сторону, противоположную движению верхней стенки.

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

При установлении течения (выход на стационарное решение) поле скорости стремится к параболическому профилю, что соответствует классическому течению Пуазейля [6]. В табл. 3 приведены значения погрешностей между скоростью, вычисленной в ходе численного эксперимента и аналитическими значениями.

2 3 5 6 8 9 11 12 14 16 -16 -12 -8-4 0 4 8 12 16 ^Ш Я u/U

а) б)

Рис. 4. Поля скорости (а) и профили скорости для различных режимов (б)

(линия - расчеты по ПК ЛОГОС, маркерами - аналитическое решение)

Таблица 3

Погрешности

Перепад давления, АР -3 -2 -1 0 1 2 3

Среднеквадратичная погрешность А 0,0014 0,0024 0,0016 1,64e-13 0,0014 0,0023 0,0018

Максимальная относительная погрешность 5,% 0,075 0,3 0,4 2,8e-10 0,176 0,15 0,08

Рисунок демонстрирует качественное и количественное согласие численного решения с точным, расхождение между аналитическими и расчетными кривыми не превышает 0,4%. Даже численный результат, в этом случае, можно считать эталонным и использовать его для сравнения с другими расчетами.

Задача 2. Ламинарное течение в канале треугольного сечения

В данной задаче рассматривается стационарное ламинарное течение несжимаемой жидкости в прямолинейной трубе с поперечным сечением в форме равностороннего треугольника (рис. 5).

Рис. 5. Сечение трубы

Геометрические характеристики расчетной области, параметры жидкости и скорость на входе выбираются таким образом, чтобы число Рейнольдса Яе=1000, для длины трубы Ь = 2 м. Сопоставление с аналитическими значениями выполняется для перепада статического давления между сечением, расположенном на расстоянии 1 м от входа в трубу и сечением, расположенном на расстоянии 0,2 м от выхода из трубы. Скорость анализируется в сечении, расположенном на расстоянии 0,2 м от выхода из трубы.

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

При стабилизированном течении в трубе с сечением в форме равностороннего треугольника перепад давления между сечениями можно вычислить по формуле [8]:

АР = 0.835 , (5)

Бг 2

где гидравлический диаметр канала определяется как Б = 2Ь/3; X - коэффициент гидравлического сопротивления.

Для ламинарного режима коэффициент гидравлического сопротивления трения X определяется по формуле[8]:

Х=64/Яе. (6)

Для стабилизированного ламинарного режима имеет место следующее распределение поперечной скорости[8]:

и ( X ) = -!- АР ( X - Ь)2Х. (7)

Рассматривались два варианта равномерных расчетных сеток: сетка 1, содержащая 150 тыс. ячеек, и сетка 2, содержащая 1200 тыс. ячеек.

В пакете программ ЛОГОС расчет выполнялся в стационарной постановке до сходимости порядка 10-6. Для аппроксимации конвективных слагаемых использовалась противопо-точная схема первого порядка точности [9].

По результатам расчета перепад статического давления на участке между двумя сечениями по длине трубы составил на сетке 1 - 2,775 Па, а на сетке 2 - 2,795 Па, что дает погрешность от аналитического значения 2,672 Па 3.8% и 4.1% соответственно.

Как и в предыдущей задаче, при установлении течения поле скорости должно стремиться к параболическому профилю. Далее на рис. 6 приведены профили скорости в сечении, расположенном на расстоянии 0,2 м от выхода из трубы, а на рис. 7 приведены изолинии скорости и поле скорости вдоль канала, полученные на сетке 2. При измельчении сетки график скорости приближается к аналитическому. Максимальная относительная погрешность наблюдается на высоте у = 0,01 м и составляет 3% для сетки 1 и 2% для сетки 2.

0.12

0.1

0.08

0.06

0.04

0.02

♦ \nalitic N сетка 1 N сетка 2

0 0 2 0 4 0 у/Н 6 0 8

Рис. 6. Профили скорости в сечении 2

Рис. 7. Изолинии скорости и поле скорости

Среднеквадратичная погрешность скорости в сетки составляет 0,000878, а для сетки 2 - 0,000413, что демонстрирует стремление среднеквадратичной погрешности к нулю. Перепад давления даже на грубой расчетной сетке не превышает 5%-ной порог погрешности, что можно считать очень хорошим результатом для применяемого численного метода и физико-математической модели.

Задача 3. Движение жидкости в квадратной каверне с движущейся стенкой

В данной задаче рассматривается квадратная каверна (рис. 8) с подвижной верхней стенкой, движущейся в своей собственной плоскости с постоянной скоростью и, боковые и нижняя стенки неподвижны. Течение считается несжимаемым, а физические свойства течения для валидации выбираются, соответствующими числу Re=1000.

На всех боковых стенках задано граничное условие «Стенка без проскальзывания», на верхней стенке - «Стенка без проскальзывания» с заданной скоростью, на передней и задней стенках - «Плоскость симметрии».

Картина течения должна содержать три вихревые структуры, что можно наблюдать и по результатам расчета (рис. 8, б): самый большой вихрь находится в центре каверны и два вихря меньших размеров в нижних углах каверны [10]. Наличие этих трех вихрей показывает, что расчетный модуль «качественно» правильно воспроизводит физику процесса.

а) б)

Рис. 8. Результаты, полученные по пакету ЛОГОС:

а - векторное поле скорости; б - изолинии

Расчеты по ПК ЛОГОС проводились в стационарной постановке до сходимости порядка 10-6 и представлены на рис. 8, а в виде векторного поля скорости.

Для расчетов использовалось два типа расчетных сеток: ортогональная равномерная сетка из четырехугольников (сетка 1: 100x100x1 = 10 000 ячеек) и нерегулярная сетка из треугольных призм (сетка 2: 22 804 ячеек). Качественно результаты расчетов по ЛОГОС совпадают с результатами из [10].

Количественные результаты расчетов данной задачи сравнивают с эталонным тестом («benchmark») из [10] по безразмерным профилям компонент скорости вдоль горизонтальной и вертикальной линий, проходящих через центр квадрата (рис. 8, б).

Рис. 9. Безразмерная горизонтальная и вертикальная компоненты скорости

В табл. 4 приведены значения среднеквадратичного отклонения расчетных значений от эталонного теста [10] для компонент скорости вдоль вертикальной и горизонтальной линий.

Таблица 4 Среднеквадратичное отклонение

Ось У X

сетка 1 2 1 2

А, % 1,4 0,007 1,33 0,013

Полученные результаты находятся в хорошем согласии с эталонными данными [10] и при измельчении сетки приближаются к ним максимально. Максимальное среднеквадратичное отклонение не превышает 1,4% для сетки 1 и 0,013% для сетки 2. Уменьшение величины погрешности на последовательности сгущающихся сеток также показывает, что численное решение сходится. Это свидетельствует о свойстве сходимости численного решения, полученного с помощью используемой физико-математической модели.

Задача 4. Безотрывное обтекание шара, задача Стокса

Наиболее старое из известных решений для ползущего течения вязкой несжимаемой жидкости принадлежит Г.Стоксу, рассмотревшему обтекание шара потоком, имеющим на бесконечности скорость их, постоянную по численному значению и направлению [11]. Течение считается несжимаемым, параметры шара и жидкости выбираются таким образом, чтобы соответствовать числу Re=0,272. Задача состоит в отыскании полной силы, действующей на шар.

Для силы сопротивления, действующей на шар при таком обтекании, справедлива формула Стокса (формула (8, а)), более же точный результат дает формула Осеена (8, б) [3]:

1 + 1 (б). (8) 16 )

Для валидации необходимо выполнить расчет в стационарной постановке до сходимости порядка 10-6. Обычно моделируется сектор с углом раствора 5 (рис. 10). Применяются следующие граничные условия: на входной границе задается постоянный расход со скоростью и, на выходной границе градиенты всех величин предполагаются равными нулю, на поверхности шара задается условие прилипания, на внешней границе расчетной области задается условие проскальзывания, на боковых стенках задано условие симметрии (рис. 10). Расчетная сетка в окрестности шара приведена на рис. 10.

Вхо,

р = (а); р = Б

Рис. 10. Граничные условия и расчетная сетка в окрестности шара

В расчете определяется полная сила, действующая на поверхность шара по формуле:

V,

Ft = F + Fp =^wAb

par

v„

+ PbAbnb'

(9)

где Ab - площадь участка поверхности стенки; vpar - составляющая вектора скорости в центре пристеночной ячейки, параллельная стенке; ^ - касательное напряжение на стенке; pb - давление на стенке; пь - внешняя единичная нормаль к стенке.

Численно полученные поле скорости и поле давления в окрестности шара представлены на рис. 11.

Поле скорости Поле давления

Рис. 11. Поле скорости и поле давления в окрестности шара

В отличие от задачи 11, в которой наблюдается отрыв потока от поверхности цилиндра в виде двух вихрей, здесь наблюдается картина стационарного безотрывного «ползущего» или медленного течения. Естественно максимальные значения давления, соответствующие минимальным значениям скоростей, наблюдаются в точке торможения потока вблизи передней части цилиндра. Такую картину течения можно наблюдать на рис. 11, на котором изображены поле скорости и поле давления в окрестности шара. Аналитическое значение силы, действующей на шар со стороны жидкости, составляет Fa = 2694,28Н, а расчетное Ff, = 2687,6Н, что составляет отличие не более чем на 0,24%.

Турбулентные течения

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

Численное моделирование турбулентности опирается на несколько подходов, основанных на описании вихревых структур различных масштабов, среди которых выделяют три основных, а именно: прямое численное моделирование (DNS, Direct Numerical Simulation) [12], моделирование крупных вихрей LES [13] (Large Eddy Simulation) и решение осреднен-ных по Рейнольдсу уравнений Навье-Стокса RANS [13] (Reynolds-AveragedNavier-Stokes).

Наиболее применимыми на практике являются RANS модели турбулентности, которые характеризуются устойчивым итерационным процессом и приемлемыми результатами для большинства практически важных типов течений. Однако они не являются универсальными и подходящими для решения широкого круга прикладных задач, и это налагает определенные требования к их применимости на практике [14]. Данного недостатка лишен «свободный от эмпиризма» DNS подход, применение которого на практике ограничено в виду потребности в огромных вычислительных ресурсах [15]. Альтернативой DNS является

«практически свободный от эмпиризма» LES подход, применение которого требует определенного качества дискретных моделей, что проявляется в существенно более мелких сетках по сравнению с RANS. Особенно это характерно для областей вблизи твердых стенок, где вихревые структуры имеют довольно малые размеры и требования к сеткам для LES приближаются к требованиям для сеток DNS [14, 15].

Понимание того, что применение LES подхода возможно только в очень отдаленной перспективе, послужило толчком для создания альтернативных моделей. Эти модели базируются на постулате о том, что в течение ближайших десятилетий на практике будут применяться RANS модели, а для моделирования крупных отрывных вихревых структур вдали от стенок, в силу намного меньших требований к вычислительным ресурсам по сравнению с DNS, можно использовать LES модели. Такие модели получили название DES или гибридных RANS-LES моделей [14, 16].

Характеризуя современное состояние полуэмпирических RANS моделей турбулентности, предназначенных для замыкания уравнений Рейнольдса, приходится констатировать, что создание универсальной модели такого типа вряд ли возможно. Однако в результате целенаправленных исследований удалось выделить две модели, имеющие наиболее высокий «рейтинг» - это модель Спаларта-Аллмараса (Spalart-Allmaras или SA модель) [17], и модель Ментера [18-19] (k - ю Shear Stress Transport модель или SST модель). Однопараметрическая модель SA исторически разрабатывалась для расчета внешних аэродинамических течений и не зарекомендовала себя при расчете задач внутренней гидродинамики, хотя модель DES, основанная на SA, дает вполне приемлемые результаты. Поэтому в дальнейшем приводимые расчеты будут основываться преимущественно на использовании SST модели. Для калибровки LES и DES моделей будет приведена отдельная задача и описаны все нюансы, которые необходимо учитывать при их валидации и калибровке.

Задача 5. Течение Хагена-Пуазейля в круглой трубе

В данной задаче рассматривается стационарное течение вязкой несжимаемой жидкости в ламинарном (Re=100), переходном (Re=3000) и турбулентном (Re=105) режимах в прямолинейной круглой трубе длиной L = 2 м, диаметром D = 0,02 м (рис. 12).

оечв-ме 1 се4»1« 2

W

2

._1_-J 0.2

L

Рис. 12. Геометрия и граничные условия

На вход в трубу подается стационарный расход жидкости, соответствующий числам Рейнольдса 100, 3000 и 105 соответственно. Для переходного и турбулентного режимов на входной границе задаются интенсивность I = 3% и длина зоны перемешивания L = 0.Ш. При стабилизированном течении в круглой трубе перепад полного давления между сечениями рассчитывается по следующей формуле [8]:

др = ^<^, (10)

О, 2

где W - средняя скорость в трубе (здесь W = w); X - коэффициент гидравлического

сопротивления трения единицы относительной длины (fj Dr=l) участка; I - расстояние между сечениями; Dr = D - гидравлический диаметр.

Для ламинарного режима коэффициент гидравлического сопротивления трения единицы относительной длины {i! Dr = l) участка определяется по закону Хагена-Пуазейля и имеет место параболическое распределение скоростей по радиусу трубы [8]:

ь = ^ (а); w (r) = ±^ I D | _ Г2

Re w 4Ц £ { 2

(б). (11)

Для переходного режима течения жидкости коэффициент гидравлического сопротивления трения определяется по данным из [8]. Для турбулентного режима коэффициент гидравлического сопротивления трения определяется по формуле Филоненко-Альтшуля, которая справедлива при Яе>4000 [8], и формуле Никурадзе, справедливой при 105 < Яе < 108 [6]:

Ь=-1-(а); 1 = 0.0032 + 0.221- Яе^237 (б). (12)

(1.8 - ^Яе-1.64)2

Профиль скорости в трубе для переходного и турбулентного режимов течения определяется как:

w (r) =u

f ff^H 'P ^

5.75' lg ^-^-+5.5

(13)

где .

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

Для сектора целесообразно использовать равномерную радиальную сетку. В расчетах по ЛОГОС использовалась расчетная сетка, содержащая 50000 ячеек (25 ячеек по радиусу трубы, 2000 ячеек по длине и 1 ячейка по углу).

Расчеты необходимо выполнять в стационарном режиме до сходимости порядка 10-6. Для аппроксимации конвективных слагаемых можно использовать любую известную схему, а в переходном и турбулентном режимах необходимо применять известные модели турбулентности. В расчетах по пакету программ ЛОГОС использовалась противопоточная схема первого порядка точности и модель турбулентности ЕЛ№ 88Т с автоматическим определением ширины пограничного слоя. В табл. 5 приведены аналитические и расчетные значения перепада давления между сечениями 1 и 2, а также значения относительной погрешности для каждого из расчетов по сравнению с аналитическими значениями.

Таблица 5

Перепад давления АР и относительная погрешность

т

Режим AP ЛОГОС, Па APэмпир, Па 5, %

Ламинарный 0,3204 0,32 0,14

Переходный 20,5 18 14,2

Турбулентный 8637,4 9230,27 6,4

8817,09 2,04

На рис. 13 приведены профили скорости в сечении 2, полученные в ламинарном, переходном и турбулентном режимах, в сравнении с кривыми, соответствующими формуле (6).

0.25 0.2 -К 0.15

\ •

\ •

0.05

0

0 0.0025 0.005 0.0075 0.01 г, м

& 3

0.0025 0.005 0.0075 0.01 г, м

Re=100

Re=3000

Re=10

5

Рис. 13. Профили скорости для различных режимов (Сплошные линии - кривые, полученные по формулам (11, б) и (13), маркеры - ПК ЛОГОС)

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

Таблица6

Среднеквадратичная погрешность скорости

Режим Ламинарный Переходный Турбулентный

А 0,000004 0,0052 0,014

Результаты, полученные по пакету ЛОГОС для ламинарного Re=100 и развитого турбулентного режимов течения Re=105, находятся в хорошем согласии с эмпирическими значениями по перепаду давления и по величине скорости. Существенное рассогласование результатов наблюдается для переходного режима. Это объясняется тем, что используемая модель турбулентности SST с автоматическим определением ширины пограничного слоя не моделирует ламинарно-турбулентный переход. Для моделирования ламинарно-турбулентного перехода необходимо использовать специализированные модели, например, SST у-Ке9 или модель, предложенную в работе [20].

Задача 6. Турбулентное обтекание пластины

В данной задаче рассматривается течение вязкого несжимаемого нетеплопроводного газа вдоль поверхности пластины. Данная задача является известным тестом, позволяющим говорить о правильности описания развития пограничного слоя с помощью различных моделей турбулентности. Для моделирования используют пластину длиной L=1 м и толщиной Ш2=0,05 мм с закруглением на носике (рис. 14).

Рис. 14. Геометрия и расчетная сетка вблизи носика

Параметры жидкости и скорость набегающего потока берутся соответствующими числу Рейнольдса = 23 • 105. Для турбулентных параметров на входной границе задаются интенсивность I = 10% и длина зоны перемешивания L=0,000504.

Для решения данной задачи целесообразно использовать блочно-структурированная расчетную сетку со сгущениями, для которой у+ ~ 1, ее фрагмент представлен на рис. 14. Валидаци-

онные расчеты необходимо проводить в стационарной постановке до сходимости порядка 10-6. Для моделирования турбулентности необходимо применять известные модели турбулентности.

Для расчетов по пакету программ ЛОГОС, использовалась расчетная сетка, число ячеек которой составляет 34 000, для аппроксимации конвективных слагаемых используется схема первого порядка точности, для моделирования турбулентности используется модель турбулентности ББТ с автоматическим определением ширины пограничного слоя (коэффициент трения дополнительно рассчитывался по моделям $Л и к-в). Для моделирования необходимо применять следующие граничные условия. На входной границе задается постоянный расход со скоростью и. На выходной границе градиенты всех величин предполагаются равными нулю. На поверхности пластины задается ГУ «Стенка без проскальзывания». На передней и задней стенке задано условие симметрии.

В расчетах определяется коэффициент трения на поверхности пластинки, который вычисляется по формуле [21]:

cf =

2т

_w

ри2

(14)

где - сила поверхностного трения.

На рис. 15, а представлено распределение вычисленного и эмпирического коэффициента трения (продольная координата нормирована на длину пластины). Вычисленный коэффициент трения на поверхности пластинки сопоставляется со значениями вычисленными по эмпирической формуле, приведенной в работе [21]:

сг = 0.0263 Яе-17. (15)

Среднеквадратичное отклонение расчетов, вычисленное по формуле (1.2), от значений, полученных по эмпирической формуле (15) составляет 0,002.

В соответствии с [21] величина безразмерной скорости и+ в пристеночной области определяется в соответствии со следующими выражениями (рис. 15, б):

- при 0 < у+ < 10, и+ = у+ ,

- при 10 < у + < 100, и += 2.441пу+ + 4.90,

+ и т

где и = —, в которой и - тангенциальная скорость жидкости и u =

и т

т

waL, xwall - касатель-Р

ное напряжение, у+ = ^ ху , в которой у - расстояние, измеряемое от стенки. На рис. 15, б

приведен эмпирический профиль скорости и профиль скорости, построенный по результатам вычислительного эксперимента.

0.01

Cf

0.008

0.006

0.004 --

0.002

о Experiment -SST --КЕ — SA

■

-- -

0.5

х, м 1

а) б)

Рис. 15. Распределение коэффициента трения и профиля скорости

На рис. 15, а демонстрируется отклонение вычисленного коэффициента трения в области «носовой» части пластинки. Это объясняется тем, что применяемые для расчета ЯЛЫ8 модели турбулентности ББТ, $Л и к-в не обеспечивают правильное описание области лами-нарно-турбулентно перехода, для моделирования которого необходимо использовать специализированные модели [20]. В остальной области пластины все используемые модели дают достаточно хорошее совпадение с эмпирическими зависимостями.

На рис. 16 представлено вычисленное поле скорости и поле давления на пластине.

а) б)

Рис. 16. Поле скорости и поле давления на пластине

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

Задача 7. Турбулентное течение в круглом канале с внезапным расширением

В данной задаче рассматривается стационарной турбулентное изотермическое течение несжимаемой жидкости в трубе со ступенчатым расширением (рис. 17) [22].

250 мм 3300 мм

*- ■*— ->

Сечение 2, г = 477.1мм i

»

1 \

50 мм

Сечение 1

Рис. 17. Круглый канал со ступенчатым расширением

На входе в расчетную область таблично задаются экспериментальные профили скорости, кинетической энергии турбулентности и скорости диссипации в соответствии с [22]. Число Рейнольдса, рассчитанное по характерной скорости на входе, составляет □ 3.2-108. Задача решается в осесимметричной постановке. Результатами численного решения являются вычисленные скорости жидкости в контрольных сечениях трубы.

В пакете программ ЛОГОС расчет проводился на двух равномерных сетках с числом ячеек ~17 тыс. и ~1,4 млн соответственно. Для моделирования турбулентного перемешивания использовалась модель ББТ с автоматическим определением ширины пограничного слоя. Для аппроксимации конвективных слагаемых применялась схема первого порядка точности.

На рис. 18 представлены безразмерные магнитуда скорости и векторное поле скоростей, полученные в ходе моделирования по ПК ЛОГОС. Особенностью рассматриваемого

течения является наличие периодической крупной вихревой структуры сразу же за ступенькой, которая впоследствии срывается вниз по потоку. На рис. 18, б демонстрируется наличие такой структуры, а также более мелкого вихря, находящегося в самом «уголке» уступа. Физико-математическая модель, подлежащая валидации, качественно должна описывать наличие обоих этих вихрей.

а) б)

Рис. 18. Магнитуда скорости и векторное поле

Количественные характеристики сравниваются с экспериментальными данными. На рис. 19 представлено сравнение расчетных и экспериментальных данных для скоростей в сечении 1 расчетные и экспериментальные данные для осевой и радиальной скоростей в сечении 3.

я. м м К-м

-Эксперимент д ЛОГОС сетка 1 х ЛОГОС сетка 2

а) б) в)

Рис. 19. Осевая (а) и радиальная (б) скорости в сечении 1 и осевая скорость в сечении 2 (в)

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

Таблица 7

Среднеквадратичное отклонение

Сечение Скорость Погрешность

Сетка №1 Сетка №2

№ 1, r = 0.05 продольная 3,22е-03 8,57е-4

№ 1, r = 0.05 поперечная 1,84е-03 1,78е-4

№ 3, r = 0.4771 продольная 1,88е-3 1,269е-3

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

Задача 8. Турбулентное течение в плоском асимметричном диффузоре

В данной задаче рассматривается стационарное турбулентное течение вязкой несжимаемой жидкости в плоском асимметричном диффузоре. Схема диффузора представлена на рис. 20, И = 0.01 м. Постановка задачи изложена в работе [23].

у

*-м-►

21-Ь 20-Ь

Рис. 20. Геометрия диффузора

На вход в трубу подается стационарный расход жидкости, соответствующий числу Рейнольдса 20 000. С целью установления развитого турбулентного профиля скорости на входе в диффузор входной участок (при х < 0, см. рис. 20) имеет длину 20И. Для решения задачи используют следующие граничные условия. На входной границе задан постоянный расход. Для турбулентных параметров на входной границе задается интенсивность турбулентности I = 3% и длина зоны перемешивания Ь = 0.1- Ь. На внешней стенке трубы задается условие прилипания. На выходной границе градиенты всех величин полагаются равными нулю. На боковых границах задается условие симметрии. Для моделирования турбулентности может использоваться любая известная модель. Расчеты должны выполняться в стационарном режиме до сходимости порядка 10-6.

В пакете программ ЛОГОС расчеты выполнялись с использованием модели турбулентности 88Т с автоматическим определением ширины пограничного слоя. Число ячеек расчетной сетки составляет 60 000. Размер пристеночной ячейки выбран, чтобы удовлетворить требованиям «высокорейнольдсовских» моделей турбулентности (30 < у+< 100). Для аппроксимации конвективных слагаемых используется схема первого порядка точности.

На рис. 21 приведены экспериментальные [23] и расчетные значения продольной составляющей скорости в различных сечениях диффузора, полученные в расчетах по пакету программ ЛОГОС. Векторное поле скорости представлено также на рис. 21. Данное течение характеризуется наличием вихревой структуры (рис. 21, б) и обширной отрывной зоны. Модуль, подлежащий валидации, должен качественно описывать физику течения. На рис. 22 для иллюстрации показано поле модуля скорости.

Юи+х

а) б)

Рис. 21. Профили скорости в плоском асимметричном диффузоре и векторное поле скорости

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

Таблица 8

Среднеквадратичная погрешность скорости

X, см 18 23 26 39 43 51 59

А 0,687 0,637 0,192 0,198 0,136 0,129 0,127

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

Задача 9. Турбулентное течение в отводах круглого и квадратного сечения В данной задаче рассматривается стационарное турбулентное течение вязкой несжимаемой жидкости в отводах круглого и квадратного поперечного сечения с углом поворота 5 = 90° (рис. 23) со следующими геометрическими характеристиками: Ьо = 1 м (диаметр для круглого сечения, ширина стороны для квадратного сечения), внутренний радиус отвода Я = 1,0 м, внешний радиус отвода Я1 = 2,0 м, длина прямолинейного участка до изгиба ¿1 = 10Ьо, длина прямолинейного участка после изгиба ¿2 = 10Ьо.

Рис. 23. Геометрическая и расчетные модели

Параметры жидкости и скорость на входе выбираются таким образом, чтобы число Рейнольдса на входе в отвод составляло 2 -105. Для решения применяются следующие граничные условия. На входной границе задается постоянный расход со скоростью Win. Для турбулентных параметров на входной границе задаются интенсивность турбулентности

I = 5% и длина зоны перемешивания Ь = 0.1- Ъ0. На выходной границе градиенты всех

величин предполагаются равными нулю. На поверхности труб задана стенка с прилипанием. Расчеты необходимо проводить в стационарной постановке до сходимости порядка 10-6. Для моделирования турбулентности можно использовать любую из моделей ЯЛЫ8 семейства.

Сравнение расчетных и экспериментальных данных проводится для величины потерь полного давления на участке между двумя сечениями, первое из которых расположено перед изгибом, а второе - в конце прямолинейного участка после изгиба, длина которого равна L2. Перепад полного давления между сечениями может быть найден по формуле [8]:

др = (а); ^ (б), (16)

где - коэффициент гидравлического сопротивления; - коэффициент местного гидравлического сопротивления; - коэффициент сопротивления трения.

Коэффициент сопротивления трения равняется сумме двух сопротивлений: = ^ , где - коэффициент сопротивления на закругленном участке; ^ - коэффициент сопротивления трения на прямолинейном участке.

Для турбулентного режима течения коэффициент потерь на трение в поворотном участке отвода может быть определен следующим образом[8]:

=0.0175 -5- X -(а); К0 = ^ = 1.5 м (б); 1=——^—— (в), (17)

0 2 ' ' (1.8 - lgR.e-1.64)2

где Ог=Ъ0 - гидравлический диаметр; X - коэффициент потерь на трение.

Коэффициент потерь на трение на прямолинейном участке отвода и коэффициент местного сопротивления определяются следующим образом[8]:

^т2 =х-Ь2 Шг (а); ^м = А: - Б! - С! (Ъ (18)

0 21

где A1 = 1- для угла поворота 5 = 90°; В1 =' , С1 = 1.

Для расчетов в пакете программ ЛОГОС использовались сетки с числом ячеек 296400 -для отвода круглого сечения и с числом ячеек 486000 - для отвода квадратного сечения. На рис. 24 показаны фрагменты этих расчетных сеток.

Рис. 24. Фрагменты расчетных сеток для отводов квадратного и круглого сечения

Расчеты по пакету программ ЛОГОС проводились в стационарной постановке до сходимости порядка 10-6. Для аппроксимации конвективных слагаемых использовалась схема первого порядка точности. Для учета турбулентного перемешивания использовалась 88Т модель с автоматическим определением ширины пограничного слоя. В табл. 9 представлены эмпирическое, полученное по формулам, и расчетные значения перепада давления. Также в таблице приведены значения У+ и относительная погрешность расчетов.

Таблица 9

Перепад давления и относительная погрешность

Отвод квадратного сечения Отвод круглого сечения

Re AP ,Па ^^ эмп' Y+ AP , Па расч 7 5, % Y+ AP , Па расч 7 5, %

2-105 5,9001 1,071 5,8263 1,25 2,153 5,8232 1,3

На рис. 25, а для иллюстрации представлено расчетное поле скорости в плоскости симметрии (слева) и фрагмент векторного поля в повороте отвода (рис. 25. б), полученные по результатам моделирования в пакете ЛОГОС.

а) б) в)

Рис. 25. Поле скорости и фрагмент векторного поля скорости для отвода круглого сечения

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

Рассчитанный по пакету программ ЛОГОС перепад давления, полученный для рассмотренного режима течения, находится в хорошем согласии со справочными данными, относительная погрешность не превышает 1.3%.

Нестационарные течения

Задача 10. Развитие ламинарного течения в трубе

В данной задаче жидкость, находящаяся в круглой прямолинейной трубе радиуса Я, до момента времени I = 0 покоится. В момент времени I = 0 внезапно возникает перепад давления йр/йх, в дальнейшем не изменяющийся по времени. Под действием сил трения и сил инерции возникает разгонное течение, которое асимптотически переходит в течение Хагена-Пуазейля с параболическим распределением скоростей. Эта задача, сводящаяся к дифференциальному уравнению Бесселя, имеет точное аналитическое решение[7]. Двумерная осесим-метричная геометрия задачи показана на рис. 26 (х - ось вращения).

Рис. 26. Течение в трубе

Для валидации расчетных моделей при решении данной задачи рассматривают нестационарное течение в трубе радиусом 1 м и длиной 100 м вязкой несжимаемой жидкости со следующими параметрами: плотность р= 1 кг/м , кинематическая вязкость у=0,01 м /с. Целесообразно рассматривать двумерную осесимметричную геометрию в виде двух моделей «Сектор» и «Труба» (рис. 27).

Модель «Сектор» Модель «Труба»

Рис. 27. Фрагменты расчетных сеток моделей «Сектор» и «Труба»

Начальные условия: по всей области трубы и = 0. Граничное условие на стенке трубы: и = 0 при всех у = Я. На левой границе задано давление р1=10Па, на правой - р2=5Па. Для задачи в постановке «Сектор» на соответствующих плоскостях задано симметричное граничное условие. Для расчета данной задачи в пакете программ ЛОГОС использовались две расчетные сетки, содержащие 2500 и 180000 ячеек соответственно. Фрагменты расчетных сеток приведены на рис. 27.

Для расчета применялись следующие разностные схемы: схема с разностями против потока для вычисления конвективных потоков и полностью неявная схема для аппроксимации по времени. Расчет проводился с постоянным шагом по времени г = 0.1.

Сравнение максимальных скоростей, реализующихся в центре трубы, для двух моделей приведены в табл. 10. Профили скоростей для различных моментов времени приведены также в табл. 10, цветные линии соответствуют расчету для модели «Труба», черные - [7]. По оси абсцисс отложено отношение поперечной координаты у к радиусу трубы Я. По оси ординат - отношение продольной компоненты скорости и к максимальной скорости ит. Здесь ит - наибольшее значение скорости, которое она достигает в середине трубы в стационарном течении Хагена-Пуазейля:

1 ип

(19)

и =———Я2

их

2 "ГТ

Роль времени на графиках играет безразмерная величина х=уг/Я . Приведены результаты для т=0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, которые соответствуют значениям времени г = 5, 10, 15, 20, 30, 40, 50 с.

Таблица 10

Значения максимальной скорости на различные моменты времени

Время г, с Т ЛОГОС , и м/с тах ' Сектор ЛОГОС итах , м/с Труба

5 0,05 0,2540 0,2513

10 0,1 0,4831 0,4783

15 0,15 0,6692 0,6647

20 0,2 0,8097 0,8058

30 0,3 0,9999 0,9971

40 0,4 1,107 1,105

50 0,5 1,167 1,166

100 1 1,239 1,239

200 2 1,244 1,243

Профили скоростей при разгонном течении в трубе (цветные линии — расчет по ЛОГОС, черные — [1])

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

Задача 11. Ламинарное нестационарное обтекание цилиндра

Данная задача является самым известным и популярным тестом для валидации программ, моделирующих ламинарное течение вязкой несжимаемой жидкости. При обтекании цилиндра потоком вязкой несжимаемой жидкости позади него образуется правильная последовательность вихрей, вращающихся попеременно вправо и влево [7]. Такая последовательность вихрей называется вихревой дорожкой Кармана. При малых числах Рейнольдса Яе < 30 обтекание тела происходит с образованием стационарной замкнутой срывной зоны в его кормовой части. При увеличении числа Рейнольдса течение за телом становится нестационарным, неустойчивым, что приводит к разрушению зоны «ползущего» течения вокруг цилиндра и отрыву вихрей поочередно то справа, то слева.

Для образования «правильной» дорожки Кармана необходимо параметры потока подобрать таким образом, чтобы получить ламинарный поток с числом Рейнольдса Re = 100-200. При значениях числа Рейнольсда ниже данного интервала дорожка может не образовываться, а при значениях выше поток станет турбулентным и вихревые структуры начнут вырождаться в вихри меньших размеров.

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

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

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

Для расчетов использовались следующие размеры цилиндра: диаметр цилиндра О = 1 • 10"4 м, высота цилиндра Н = 5 • 10"6 м, лобовая площадь: £ = Н • О = 5 • 10"10 м2.

Рис. 28. Геометрия и сетка расчётной области

Расчеты проводились с первым и вторым порядком аппроксимации по времени. Поле скорости и давления представлены на рис. 29. В отличии от задачи с ползущим характером течения, здесь наблюдается поочередный срыв вихревых структур с поверхности цилиндра, образующих вниз по потоку «вихревую» дорожку Кармана. Вихри имеют периодичную

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

мдв и\м/ * р

0.00 0 06 0.13 0.19 0.25 0.32 0.38 0.44 0.51 0.57 0.63 -0.15 -0.12 -0.09 -0.05 -0.02 0.01 0.04 0.08 0.11

Рис. 29. Поля распределения скорости давления

Количественные характеристики оцениваются с помощью расчета относительно1 погрешности числа Струхаля, коэффициента лобового сопротивления и периода колебаний гидродинамических величин [7]. Результаты валидации пакета программ ЛОГОС для данной задачи представлены в табл. 11.

Таблица 11

Результаты расчета и относительная погрешность

Параметр Эксперимент ЛОГОС 8, %

Число Струхаля 0,1895 1-й порядок 0,1825 3.68

2-й порядок 0,1865 1.57

Коэффициент лобового сопротивления 1,4 1-й порядок 1,4999 7.14

2-й порядок 1,4748 5.34

Период колебаний подъемной силы 0,00124 1-й порядок 0,00129 4.03

2-й порядок 0,00126 1.61

Попадание результатов в 5%-ный интервал погрешности можно считать удовлетворительным результатом. Именно допустимый интервал погрешности в 5% считается нормальным при количественной оценке результатов расчета. Выход численного решения из этого интервала означает недостаточную точность используемой физико-математической модели или численной схемы. Точность модели можно повысить путем измельчения расчетной сетки или повышением порядка численной схемы (хотя и это может не сработать и нужно более детально анализировать задачу и применяемую математическую модель).

02

02

Рис. 30. График зависимости скорости V от времени в точке слежения

Количественные оценки колебательного режима оцениваются по колебаниям величин

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

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

Задача 12. Поперечное обтекание трубного пучка

В данной задаче рассматривается нестационарное течение воды через канал с пучком труб в трех различных конфигурациях. Для этой задачи можно провести аналогию с предыдущей - здесь моделируется обтекание системы из нескольких цилиндров в различных конфигурациях. Геометрия конфигураций для трех таких систем представлена на рис. 31 [24].

а) б) в)

Рис. 31. Расположение цилиндров в трех различных конфигурациях:

а - расположение в линию; б - шахматный порядок; в - шахматный порядок

Размеры расчетной области в направлении течения выбираются таким образом, чтобы от начала области до центра первого стержня расстояние составляло 10 ё, а от центра последнего до конца области - 12,5 ё (здесь ё - диаметр цилиндра).

Для валидации необходимо решить нестационарную задачу. Жидкость рассматривается как несжимаемая со скоростью потока на входе 0,93 м/с. Плотность воды берется равной 997,561 кг/м3, молекулярная вязкость 8,8871е-4 кг/(м-с), что соответствует Яе « 104.

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

Размерность сетки для различных случаев расположения:

1) 716977 сеточных элементов для расположения в линию, параметры решетки 3,6х2,1 см (рис. 31, а);

2) 591785 сеточных элементов для расположения в шахматном порядке, параметры решетки 3,6х2,1 см (рис. 31, б);

3) 470018 сеточных элементов для расположения в шахматном порядке, параметры решетки 3,6х1,6 см (рис. 31, в).

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

гаемых используется схема первого порядка точности, шаг по времени выбирается равным Д^=3,240-5а Для случая расположения пучка стержней в линию и для расположения в шахматном порядке с параметрами решетки 3,6х2,1 см конечный момент времени расчета ¿=1 с. Для случая расположения в шахматном порядке с параметрами решетки 3,6х1,6 см конечный момент времени расчета ¿=0,4 с, что связано с меньшим размером расчетной области.

Граничные условия задаются в соответствии со схемой, представленной на рис. 32.

Рис. 32. Фрагмент расчетной сетки для расположения в шахматном порядке стержней, параметры решетки 3,6x2,1 см

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

Рис. 33. Сечения расчетных областей

На рис. 33-35 выборочно представлены зависимости продольной компоненты средней скорости для различных сечений для трех вариантов расположения стержней. Численное решение сравнивается с экспериментальными данными[24].

Рис. 33. График продольной составляющей средней скорости для расположения в линию (маркеры - эксперимент, линия - результат по ПК ЛОГОС)

Рис. 34. График продольной средней скорости для шахматного порядка 3,6x1,6 см (маркеры - эксперимент, линия - результат по ПК ЛОГОС)

Рис. 35. График продольной средней скорости для шахматного порядка 3,6х2,1 см (маркеры - эксперимент, линия - результат по ПК ЛОГОС)

В табл. 12 приведены значения среднеквадратичного отклонения А продольной компоненты средней скорости для сечений х/ё=0, 2.95, 5.9, 8.4. На некоторых линиях наблюдается некоторое отклонение от эксперимента, что может быть связано с недостаточным разрешением сеточной модели или неточностью описания переходных областей течения.

Таблица 12

Значения среднеквадратичного отклонения

х / d 0,0 2.95 5.9 8.4

в линию

А 0,020604 0,02208 0,014571 0,024826

шахматный порядок 3,6х1,6 см

А 0,0107 0,2119 0,0469 0,0714

в шахматный порядок 3.6х2.1 см

А 0,015798 0.054710 0,062805 0,041995

На рис. 36-38 показаны сцены мгновенного распределения скоростей и турбулентной кинетической энергии в расчетных областях.

Рис. 36. Распределение скорости и турбулентной кинетической энергии на момент времени ^ = 1 с (в линию)

Рис. 37. Распределение скорости и турбулентной кинетической энергии на момент времени ^ = 0,4 с (шахматный порядок 3,6х1,6 см)

Рис. 38. Распределение скорости и турбулентной кинетической энергии на момент времени ^ = 1 с (шахматный порядок 3,6х2,1 см)

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

Задача 13. Вырождение однородной изотропной турбулентности

Как отмечалось ранее, одним из основных подходов к моделированию нестационарной турбулентности в настоящее время является использование DES модели. Для успешного использования моделей такого класса важны свойства выбранной схемы дискретизации, от которой напрямую зависит качество DES моделирования: схема должна иметь малую диссипацию и обеспечивать устойчивость счета. Выбор оптимальной схемы дискретизации для конвективных потоков - одна из основных проблем при моделировании течений с помощью моделей DES. Схема должна, с одной стороны, иметь малую диссипацию, т.е. порождать как можно меньше численной диффузии, с другой стороны, она должна обеспечивать стабильный счет на произвольной неструктурированной сетке. Требование малой численной диффузии проистекает из того факта, что при DES моделировании в основной области, где присутствуют вихревые структуры, используется модель LES, и для правильного описания эволюции данных структур необходимо, чтобы численная диффузия схемы была меньше, чем диффузия, которая обеспечивается подсеточной вязкостью модели LES. Поэтому при вали-дации методов DES особое внимание должно уделяться схемам дискретизации. В настоящей задаче проводится анализ применимости для DES моделирования следующих схем:

• противопоточная схема с линейной интерполяцией (LUD);

• схема QUICK;

• центрально-разностная схема (CD);

• центрально-разностная схема с линейной реконструкцией (CDW);

• схемы семейства NVD;

• гибридные схемы (перечисленные ранее схемы, смешанные с противопоточной схемой для увеличения монотонности).

Схема LUD - противопоточная схема с использованием реконструкции величины на грань с линейной интерполяцией [25, 26]. Схемы CD и CDW- наименее диссипативные схемы, однако они неустойчивы [25, 27], и при наличии больших градиентов их использование приводит к осцилляциям в поле решения. Схема QUICK [27] является схемой смешанного типа, она сочетает в себе противопоточную схему первого порядка совместно с градиентной

поправкой и схему CD. Такой подход позволяет убирать диссипативные ошибки первого порядка и подавлять осцилляции второго порядка. Схемы семейства NVD основаны на использовании диаграмм нормализованной переменной и критерия ограниченности CBC [28]. В своей изначальной форме построение функции нормализованной переменной основано на использовании структурированной сетки. Наиболее удобным обобщением схем NVD на произвольные неструктурированные сетки рассмотрено в работе [29], там же предложена сама схема семейства NVD - схема GAMMA. Схема GAMMA имеет в основе своей схему CD, но использует функцию нормализованной переменной для применения критерия ограниченности CBC. Схема GAMMA более диссипативна, чем CD, однако приводит к ограниченному решению искомой величины. Также схема имеет один свободный параметр, который определяет вклад противопоточной составляющей в результурующее значение величины на грани и, в конечном счете, влияет на быстроту сходимости.

Гибридные схемы представляют собой линейную комбинацию схемы высокого и низкого порядков [26], что приводит к увеличению монотонности решения. Гибридную схему, например, для схемы CD и UD, можно записать следующим образом:

<k,Sibr = y<t>F,CD + (1 - У)ФеДВ , (20)

где у - коэффициент смешения. В настоящей работе рассматриваются гибридная схема на основе CD и UD, и схема на основе CDWи LUD, с постоянным коэффициентом смешения у = 0,9.

Оценка диссипативности схем проводится путем решения задачи о вырождении однородной изотропной турбулентности [30]. Однородная изотропная турбулентность - достаточно редко реализующееся в природе явление, но его моделирование помогает с хорошей точностью узнать насколько диссипативна та или иная схема дискретизации конвективных потоков. По диссипативности можно оценить пригодность каждой схемы для DES расчетов в области LES.

Моделирование задачи проводится в расчетной области размером 2л>2л>2л: с количеством ячеек 64-64-64 по всем трем направлениям соответственно. По всем трем пространственным координатам используется периодическое граничное условие. По заданному начальному энергетическому спектру формируется случайное поле скоростей, которое используется в качестве начального условия. Далее производится нестационарный расчет задачи до момента времени 0,87 с, после чего результаты сохраняются и переводятся обратно в энергетический спектр, который сравнивается с экспериментальными данными [30]. По поведению кривой энергетического спектра можно судить о качестве каждой схемы - занижение энергии в том или ином интервале говорит о ее диссипативности.

Для сравнения схем по уровню численной диссипации в расчете использовалась модель турбулентности LES со значениями констант Смагоринского Cs = (0;0,1;0,2). Для дискретизации по времени применяется схема Кранка-Николсона, шаг по времени соответствует числу Куранта равному единице. На каждом шаге по времени проводится 10 внутренних итераций решателя, что обеспечивает сходимость невязок основных физических величин до S =10-8 по норме L1.

Далее представлены результаты расчета задачи для схем LUD и QUICK. Данные схемы достаточно диссипативны и в области высоких частот при любых значениях константы наблюдается занижение энергии (рис. 39). Даже обнуление константы Cs, а, следовательно, и подсеточной вязкости, не приводит к правильной картине.

На рис. 40 представлены результаты для схемы GAMMA. Расчет проведен для рекомендованного максимального значения параметра схемы у=0,1 [25]. При таком значении у доля противоточности в схеме минимальна и она наименее диссипативна. Также проведены расчеты для схемы у=0,2, которая более применима для неструктурированных сеток в областях со сложной геометрией. По полученным результатам можно сделать вывод, что схема GAMMA, так же, как и QUICK дает занижение энергии в области высоких частот при любых значениях константы Cs.

Рис. 39. Энергетические спектры, схема LUD и QUICK

1.0Е-01

1.0Е-02

1.0Е-03

1.0Е-04

Е(к) Схема GA ММА -0.2

\

1.0Е+00

1.0Е+01

1.0Е+02

Эксп.-Cs = 0.0 —Cs = 0.1 —Cs = 0.2

Рис. 40. Энергетические спектры, схема GAMMA

Расчеты показывают, что наименее диссипативна схема CD, которая при Cs = 0,2 правильно описывает эволюцию вихревых структур (рис. 41). Однако устойчивый счет при использовании данной схемы на неструктурированных сетках возможен только при использовании гибридной схемы с коэффициентом смешения у=0,9 и меньше (формально обозначим схему при у=0,9: 0,9CD+0.1UD). Далее представлены результаты расчета для 0.9CD+0.1UD. По результатам видно, что добавление 10% доли противопоточности уже заметно повышает диссипативность схемы, однако для нее всё-таки удается подобрать оптимальное значение константы Cs = 0,1.

1.0Е-01

1.0Е-02

1.0Е-03

1.0Е-04

Е(к) Схема 0 .9CD+0.1UD

/ >

\

1.0Е+00 1.0Е+01 1.0Е+02 ■ Эксп.-Cs = 0.0-Cs = 0.1-Cs = 0.2

Рис. 41. Энергетические спектры, схема CD и 0,9CD+0,1LD

Для иллюстрации на рис. 42 приведены поля модуля скорости, давления и турбулентной вязкости на конечный момент времени - 0,87 с, полученные в результате расчета для 0,9СВ+0,1ЦО. Эволюция вихревых структур различных масштабов показана на рис. 43.

а)

б) в) Рис. 42. Поля модуля скорости (а), давления (б) и турбулентной вязкости (в) на 0,87 с

0,27 с

0,58 с

Рис. 43. Эволюция турбулентности

0,87 с

Таким образом, наименее диссипативными вариантами схем являются CD и 0,9CD+0,1UD: при DES расчете в области LES они смогут обеспечить правильное описание эволюции вихрей и передачу энергии от крупных вихрей к мелким. Схемы QUICK и GAMMA более диссипативны, чем CD и 0,9CD+0,1UD. Для метода конечных объемов величина численной диссипации зависит также и от типа элементов, из которых составлена расчетная сетка. Однако, на сетках различного типа тенденция давать численную диссипацию у схем будет такая же: на одной и той же сетке LUD будет диссипативней, чем QUICK, а QUICK диссипативней, чем CD и 0,9CD+0,1UD. Поэтому для дальнейших расчетов с использованием DES модели будем использовать только схемы CD и 0,9CD+0,1UD, предварительно проведя для них калибровку констант.

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

а) б)

Рис. 44. Энергетические спектры, схема CD (а) и 0,9CD+0,1UD (б)

Для схемы CD оптимальные значения констант составляют £¿£-=0,61, Сю=0,78 для схемы 0,9CD+0,1UD, значения констант должны быть примерно в два раза ниже: £¿¿=0,305, Ою=0,39. Именно эти значения констант для каждой схемы и должны использоваться в расчетах по пакету программ ЛОГОС.

Задача 14. Развитое турбулентное течение в канале

Данная задача является хорошим тестом для валидации вихреразрешающих моделей типа DES. В данной задаче рассматривается развитое течение несжимаемой жидкости в плоском канале с высотой 1 м при числе Рейнольдса, построенном по высоте канала и динамической скорости, равном ReT = 800, что соответствует условиям, при которых в работе [12] выполнен DNS расчет данного течения. Для определения осредненных характеристик течения вначале необходимо выполнить нестационарный расчет до выхода решения на стати-стически-установившийся режим. Далее продолжать расчет с вычислением осредненной горизонтальной скорости и осредненных компонент тензора турбулентных пульсаций.

В качестве модели турбулентности необходимо использовать либо LES, либо DES модели. В расчетах по пакету программ ЛОГОС, представленных далее, использовались три модели DES (DES-ks, DES-SA, DES-SST) с различными схемами дискретизации конвективных слагаемых в уравнениях переноса импульса и соответствующими наборами констант из табл. 13.

Таблица 13

Оптимальные значения эмпирических констант

CD CDW 0,9CD+0,1UD 0,9CDW+0,1LUD

Cks 1,20 1,18 0,61 0,60

CSA 0,62 0,60 0,33 0,31

CSST-km 0,78 0,72 0,39 0,38

CSST-ks 0,61 0,58 0,31 0,30

Для расчета использовалась сетка, содержащая сгущения к областям твердых поверхностей. Общее количество ячеек составляет 1,2 млн. На твердых стенках компоненты скорости, турбулентная вязкость и кинетическая энергия турбулентности полагаются равными нулю. На остальных границах для всех переменных используется условие периодичности. Для того чтобы обеспечить движение, в уравнение переноса продольной составляющей импульса вводится дополнительный член, равный градиенту давления при установившемся течении в канале: dp / dx = 1.721 Па/м.

На рисунке в качестве иллюстрации представлены мгновенные поля скоростей, полученные в результате моделирования по трем разным подходам к моделированию турбулентности (RANS, DES и LES).

RANS LES DES

Рис. 45. Мгновенное поле скорости

На рис. 46-47 представлены профили средней скорости и разрешенных компонент тензора напряжений Рейнольдса для каждой схемы в сравнении с результатами, полученными при DNS расчете [12].

DES КЕ

и+

^Г -Г)МЧ Mozer)

-CDW 0.9СС -CD -0.9СС W+0.1LUD +0.1UD

V+

DES SA

-DNS (Mozer]

CDW CD —0.9CD+0.1UD 0.9CDW+0.1LUD

у+

1

0,8 0,6 0,4 0,2 U'VO -0,2 -0,4 -0,6 -0,8 -1

1 0,8 0,6 0,4 0,2 u'vO -0,2 -0,4 -0,6 -0,8 -1

DES KE

'm

I 0 2 0 4 0 6 0 8

□NS (Maze r et all)

0,9CD+0,1 UD

0.9CDW+0 1CD

CDW

DES SA

J

I 0 2 0 4 0 6 0 8

......DNS (L lozer)

CD

-CDW

-0.9CDV V+0,1CD

Рис. 46. Профиль скорости и напряжений Рейнольдса (DES-ks, DES-SA)

Рис. 47. Профиль скорости и профиль напряжений Рейнольдса (DES-SST)

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

Смешанная и естественная конвекция

Задача 15. Смешанная конвекция в вертикальном плоском канале

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

В данной задаче рассматривается двумерное стационарное ламинарное течение в длинном вертикальном канале (рис. 48). Температуры стенок канала равны Т^ и Тс (Тъ > Тс). На входе в канал задается однородный поток, температура которого равна полусумме температур стенок. На выходе из канала задается постоянное давление и ставятся «мягкие» условия (д*/ду=0) для скорости и температуры. Ускорение силы тяжести g направлено в сторону, противоположную оси у.

Тс

Т(1

ТТТТ

Рис. 48. Схема вертикального канала

В качестве масштаба длины принята ширина канала Ь, характерный перепад температуры равен АТ = Ть - Тс, характерной температурой считается полусумма температур стенок канала Т = (Тн+Тс)/2, характерной скоростью - среднерасходная (она же входная) скорость V. Безразмерные граничные условия формулируются следующим образом: ы=у=0, 9 =- 0.5 при х=0; ы=у=0, 9=+ 0.5 при х=1; V = 1, 9 = 0 при у=0 (вход) р =0, дu/дy=дv/дy=д9/дy=0 при у=Н (выход).

Размеры области, параметры среды и граничные величины в расчетах приняты такими, чтобы число Рейнольдса Яе=100, число Прандтля Рг=1 и число Грасгофа 0г=0;3600;7200.

На участке установившегося течения, где линии тока параллельны оси канала, профили поперечной составляющей скорости v(x) и избыточной температуры 9(х) не зависят от продольной координаты у и описываются следующими уравнениями [3]:

Ц = Яе ^ - ^ .9, йх2 йу Яе

*,

й 29

йх 2

= 0.

(21)

Продольный градиент давления ёр /ёу есть неизвестная константа, для определения которой служит условие нормировки скорости (сохранения расхода):

0.5

уёх = 1.

(22)

-0.5

С учетом вышеприведенных граничных условий на стенках канала решение данной задачи можно записать в следующем виде [3]:

9 = х-1, V = ** (х - х3)+Г 6* (х - х2 ) ^ = -™, * э(23) 2 6 V ' \ 4 1 V р ду Яе Яе

При значениях параметра плавучести В=72 производная (^/ёх)х=0 обращается в ноль; при В > 72 вблизи левой стенки канала (х=0) появляется область нисходящего течения V < 0). Отметим, что в рассматриваемом случае величина градиента давления ёр /ёу не зависит от наличия сил плавучести и соответствует известному закону сопротивления для чисто напорного течения в плоском канале (течение Пуазейля; коэффициент сопротивления ^=24/Яе).

Расчеты в пакете программ ЛОГОС проводились на равномерной сетке размером 50x3000^1 ячеек. Для решения задачи задаются следующие граничные условия: на вертикальных стенках - стенка без проскальзывания с заданной температурой, на нижней границе - вход, на верхней - выход, на передней и задней поверхностях расчетной области - плоскости симметрии.

Задача в пакете программ ЛОГОС решается с использованием модели вязкой несжимаемой жидкости. Для учета сил плавучести (термосжимаемости) используется модель Буссинеска. [3]. Поля температуры и скорости потока, полученные при числе Грасгофа Ог = 7200 (что отвечает «критическому» значению параметра плавучести В = 72), представлены на рис. 49.

Рис. 49. Поле скорости и температуры при числе Грасгофа Сг = 7200

На рис. 50 представлено сравнение профилей скорости, полученных при различных значениях параметра плавучести В=Ог/Яе, с решением (23) для установившегося течения. На рис. 51 представлены графики изменения градиента давления ёр /ёу вдоль канала. Представлено точное решение (12/Яе = 0.12) и решение, полученное по ПК ЛОГОС.

1.8 1.6 1.4 1.2 1 1 >0.8 0.6 0.4 0.2 0

в N

.........-1 f // ........... д

0.2

0.4

0.6

0.8 1 х, м

Рис. 50. Профили скорости V Рис. 51. Изменение градиента давления

в вертикальном канале вдоль канала

Аналитическое решение - маркеры, ЛОГОС - линии

Как видно из рис. 50, соответствие расчетного и аналитического решений практически полное. Величина относительной погрешности 5<0,1%, а величина среднеквадратичного отклонения А =0,00043. Представленные на рис. 51 графики изменения градиента давления йр /йу вдоль канала свидетельствуют, что длина начального участка быстро увеличивается с ростом параметра плавучести. При этом, если для чисто напорного течения (В=0) градиент давления перестал меняться при у~10, то для случая В=72 установившийся режим течения, строго говоря, не был достигнут до выхода из канала. Тем не менее, для всех рассмотренных режимов правильная асимптотика йр /йу при удалении от входа в канал очевидна.

Задача 16. Ламинарная естественная конвекция в квадратной полости

В данной задаче рассматривается стационарное ламинарное течение вязкой несжимаемой жидкости в квадратной полости с длиной стороны Ь, находящейся в поле действия силы тяжести. Температура правой и левой вертикальных стенок соответственно равна Ть и Тс (Ть > Тс); горизонтальные стенки теплоизолированы (рис. 52) [31, 32]. Для решения задачи применяются следующие граничные условия. Вертикальные стенки моделируются как стенки с прилипанием с заданными температурами Тс и Ть. Горизонтальные стенки рассматриваются как теплоизолированные стенки с прилипанием. На внешние плоскости, повторяющие исходную двумерную геометрию, накладываются граничные условия «плоскость симметрии».

0.6

'0.4

0

Т 1 р\ —1 I

тг1 г< / > ^ 06

-80 -40

40

80. м/с

Рис. 52. Геометрия

Рис. 53. Горизонтальная составляющая скорости на средней линии Сплошная линия - ЛОГОС, маркеры - [32]

В качестве масштаба длины принята сторона квадрата Ь, характерный перепад температуры равен АТ = Ть - Тс, характерной температурой считается температура холодной стенки Тс. Рассматриваются две постановки задачи, в которых параметры среды и граничные величины заданы таким образом, чтобы число Релея было равно Яа = 104 и Яа = 106 соответственно. Число Прандтля Рг=0,71.

Расчеты по пакету программ ЛОГОС выполнялись в стационарной постановке до сходимости порядка 10-6. Для аппроксимации конвективных слагаемых используется схема первого порядка точности. Естественная конвекция моделируется путем учета влияния массовых сил с помощью модели Буссинеска. В расчетах используется равномерная сетка, содержащая 80^80x1 ячеек, такая же, как и в работе [32].

Для сравнения с результатами работы [32], полученные в расчетах скорости приводятся к масштабу а/Ь, где а=Х/рСр - коэффициент температуропроводности. На рис. 52 представлен профиль горизонтальной составляющей скорости на средней линии (х=Ь/2) для сетки 80x80. Среднеквадратичное отклонение расчета по ПК ЛОГОС от результатов в работе [32] составляю,т 0,29 для Яа = 104 и 0,56 для Яа = 106.

В табл. 14 приведены результаты, полученные в ПК ЛОГОС, в сравнении с результатами, полученными различными авторами и заимствованными из работы [31]. В таблицах

величины итах, Утах - значения максимальных скоростей на линиях х = Ь/2, у = Ь/2 соответственно; № = дЬ/ХАТ - число Нуссельта, где q - средний тепловой поток.

Таблица 14

Результаты расчета

U nax V r max N u Источник

Ra=104 Ra=106 Ra=104 Ra=106 Ra=104 Ra=10 6

16,1798 64,6912 19,6177 220,8331 - - Mayne [32]

16,1 65,4 19,9 228 2,084 8,743 Manzari [32]

16,2 65,33 19,51 216,75 2,243 8,8 De Vahl Davies [32]

16,149 71,532 19,427 223,589 2,25 9,096 ЛОГОС

На рис. 54-55 представлены векторные поля скоростей и поля температуры, полученные в расчетах по пакету ЛОГОС на сетке 80x80 ячеек.

Ra=104 Ra=106

а) б)

Рис. 54. Векторное поле скорости (сетка 80х80)

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

Ra=104 Ra=106

а) б)

Рис. 55. Поле температуры (сетка 80х80)

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

Задача 17. Развитая турбулентная конвекция

Данная задача похожа на предыдущую. Здесь исследуется развитое турбулентное конвективное течение воды в подогреваемой снизу кубической полости с длиной ребра 250мм [33]. Вертикальные стенки куба теплоизолированы плитами из пеноплекса, горизонтальные выступают в качестве теплообменников. Средняя температура воды в полости 250С, что соответствует среднему числу Прандтля Pr = 6,1. Заданное значение перепада температуры между теплообменниками ДТ = 20 0С (нижний теплообменник имеет температуру Т1 = 350С, верхний - Т2 = 150С), что соответствует числу Релея Ra = 6,1x10 . Результаты расчетов сравнивались с экспериментальными данными, полученными в [33].

Расчеты по пакету программ ЛОГОС выполнялись в нестационарной постановке. Расчет проводился в два этапа: расчет развития турбулентной конвекции, за которым следовал расчет развитой турбулентной конвекции с осреднением расчетных величин. Для аппроксимации конвективных слагаемых используется схема второго порядка точности - CD, для дискретизации по времени использовалась схема второго порядка Адамса-Бешфорта. В расчетах используется равномерная сетка, содержащая 150*150x150= 3375000 ячеек. На рис. 56 представлены средние поля скорости, полученные в эксперименте и в расчетах.

в)

г)

Рис. 56. Средние поля скорости:

эксперимент: а - модуль скорости; б - векторное поле средней скорости; расчет Логос: в - модуль скорости; г - векторное поле средней скорости

В эксперименте и в расчетах структура средних полей скорости очень близка. Хорошо видно, что в полости возникает крупномасштабный вихрь, соизмеримый с размерами полости. Помимо основного вихря на средних полях существуют вторичные вихри меньшего масштаба, локализованные в противоположных углах кубической полости. Для количественного сравнения были построены проходящие через центр полости при фиксированной координате (125 мм) профили горизонтальной и вертикальной компонент скорости (рис. 57). Оба профиля демонстрируют хорошее совпадение, повторяя не только структуру потока, но и воспроизводя значение скорости в пограничных слоях.

а)

б)

Рис. 57. Профили средней по времени скорости (синяя линия - расчет; красная штриховая линия - эксперимент):

а - профиль вертикальной компоненты скорости); б - профиль горизонтальной компоненты скорости

Для горизонтальной компоненты скорости Среднеквадратичное по профилю отклонение расчета от эксперимента составило 0,17 мм/с. Для вертикальной компоненты скорости аналогичные оценки дают среднеквадратичное отклонение в 0,62 мм/с.

Течения с теплообменом

Задача 18. Течение в канале с обратным подогреваемым уступом

Данный тест является одним из сложнейших в вычислительной гидродинамике. В нем рассматривается течение вязкого несжимаемого газа в канале с расширением в виде обратного уступа и подогреваемой за ним стенкой. Для данной задачи имеются экспериментальные данные, представленные в [34]. Схема геометрии расчётной области представлена на рис. 58. Высота уступа Н=0,038 м. На рис. 59 показана схема течения потока за обратным уступом.

Рис. 58. Схема геометрии расчётной области

Рис. 59. Схема течения потока за обратным уступом

Физика течения содержит обширные отрывные зоны, зону циркуляции потока, состоящую из двух вихрей различных масштабов, зону присоединения большого вихря, зону отрыва большого вихря и другие особенности. Для многих моделей турбулентности данный тест остался недосягаемым. Из RANS моделей турбулентности наиболее адекватно это течение описывает модель SST, но все равно имеет проблемы с описанием коэффициента трения. Для моделирования турбулентности необходимо использовать модели, основанные на вих-реразрещающих подходах [9], которые дают более точные результаты.

Для моделирования параметры среды и граничные величины выбираются так, чтобы Re=27867 и Pr=0.7. Применяются следующие граничные условия. На входной границе задается постоянный расход со скоростью U. На выходной границе градиенты всех величин предполагаются равными нулю. На внешней границе расчетной области задается условие проскальзывания. На нижней стенке за уступом задан тепловой поток. На боковых стенках задано условие симметрии.

По результатам расчета определяются коэффициент трения Cf на нижней подогреваемой стенке (канала с расширением в виде уступа) и коэффициент Стантона St.

Коэффициент трения и число Стантона вычисляются по формулам:

С, = % »(х) = ,

ри и-р-Ср

где т№ - сила поверхностного трения; а/х\ _ Ч - локальный коэффициент теплоотда-

( ) Т (х) - Т

чи; Tb - температура потока; Cp - теплоемкость; q - тепловой поток на подогреваемой стенке; ^ - расчетная температура на стенке.

Для расчета целесообразно использовать блочно-структурированную сетку со сгущениями (на входе в зону расширения, в области присоединения потока), для которой У+~1 (рис. 60). При моделировании по пакету программ ЛОГОС число ячеек сетки составило 18930.

Рис. 60. Сеточная модель

В расчетах по пакету программ ЛОГОС используется модель SST с автоматическим определением ширины пограничного слоя, а также DES модель. Расчет выполнялся в стационарной постановке до сходимости порядка 10-6. Для аппроксимации конвективных слагаемых используется схема первого порядка точности.

На рис. 61 для иллюстрации показано поле распределения скорости, полученное в расчете по пакету ЛОГОС.

MAG_UVW *

2.11096-005 2.9857 5.9714 8.9571 11.9428

Рис. 61. Поле распределения скорости

Расчет показал формирование за обратным уступом зоны торможения потока, отмечено образование большого 1 и малого 2 вихрей (рис. 62).

Рис. 62. Течение за обратным уступом

Как и в задаче 7, большой вихрь, в последствие, срывается вниз по потоку. Именно зона «нестационарного» отрыва и является очень проблематичной для моделирования с применением RANS моделей. На рис. 63 представлено распределение коэффициентов Cf и St

г _ Х ~ Хтш

расчетное и экспериментальное от нормированного значения хпагш .

Хтах хтт

х norm

а) б)

Рис. 63. Распределение коэффициента трения (а) и числа Стантона (б) в сравнении с [34] (сплошная линия расчет по ЛОГОС)

Распределение коэффициента трения (рис. 63, а) показывает существенно заниженные значения вблизи зоны отрыва и эти вычисленные значения с помощью наиболее удачной RANS модели турбулентности. Исправить ситуацию позволяет применение гибридной DES SST модели. На рис. 64 представлены результаты численного решения данной задачи с помощью IDDES модели и сравнение результатов с SST моделью [9].

Схема CD

0.003

0.002

<0.002

■ Эксп. -IDDES -SST

Схема 0.9CD+0.1UD

0.003 -т 0.002

0.002

■ Эксп. -IDDES -SST

Рис. 64. Сравнение коэффициента трения по моделям IDDES и SST

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

Для представленных ранее величин было рассчитано среднеквадратичное отклонение. Для коэффициента трения среднеквадратичное отклонение численного решения от эксперимента составляет 2.1 -10"5, для числа Стантона 7.8 -10"5. На основании проведенного анализа полученных результатов, по распределениям коэффициента трения и числа Стантона, можно сделать вывод, что результаты расчета качественно и количественно согласуются с экспериментальными данными (отклонение точки присоединения не превышает 5%). Аналогичные значения должны получаться и для вновь разработанных программ, подлежащих валидации.

Задача 19. Охлаждение твердотельных блоков в плоском канале

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

В данном тесте рассматривается задача охлаждения шести твердотельных блоков с заданным объемным тепловыделением, установленных в плоском прямоугольном канале, ламинарным потоком воздуха при вынужденной конвекции [35]. Геометрия расчетной области схематично представлена на рис. 65.

У

н ^т' шт мш --ь.

5.6Н НИ х <-►

1-=20Н

Рис. 65. Схема геометрии расчётной области

Плоский прямоугольный канал имеет высоту Н = 1 м и длину Ь = 20Н. Внутри него, вплотную к нижней стенке, расположены шесть одинаковых твердотельных блоков длинной Н и высотой И = 0.25Н . Левая грань первого блока находится на расстоянии 5.6Н от входа. Промежутки между блоками одинаковы и равны Н. Каждый блок равномерно нагревается объемным источником тепла мощностью Q. Стенки канала неподвижны, непроницаемы и теплоизолированы. Эффекты плавучести и вязкой диссипации не учитываются. Параметры жидкости выбираются так, чтобы число Рейнольдса Яе=100, а число Прандтля Рг=0.7.

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

В пакете программ ЛОГОС расчет выполнялся в стационарной постановке на расчетной сетке из 4000 ячеек. Для аппроксимации конвективных слагаемых использовалась схема второго порядка точности. Размеры сетки выбираются так, чтобы границы твердотельных блоков точно совпадали с гранями ячеек.

В табл. 15 представлены параметры жидкости и твердого тела.

Таблица 15

Параметры жидкости и твердого тела

Параметр Значение Размерность

Параметры жидкости

Плотность жидкости, Pf 1 кг/м3

Молекулярная вязкость ц 0,01 кг/(м-с)

Коэффициент теплопроводности X 0,0142857 Вт/(мК)

Теплоемкость при постоянном давлении Ср 1 Дж/(кгК)

Скорость набегающего потока 1 м/с

Параметры твердого тела

Плотность твердого тела, р3 100 кг/м3

Объемный источник тепла Q 71,142857 Вт/(м3)

Коэффициент теплопроводности X 71,142857 Вт/(мК)

Теплоемкость при постоянном давлении Ср 1 Дж/(кгК)

На рис. 66 приведено поле температуры, полученное по ПК ЛОГОС. Видно, что твердотельные блоки разогреваются за счёт объемного источника тепла. Расположенные ближе к входу блоки имеют меньшую температуру, поскольку обтекаются жидкостью с меньшей температурой. На рис. 67 приведены профили температур вдоль прямой У = 0.1м, полученные по ПК ЛОГОС и в работе [35].

Рис. 66. Поле температур

40 30

Ьо

н

10

♦ Bencafada I

--ЛОГОС ' МП

nfi

г Ур

► i

0 5 10 15 20

R, м

Рис. 67. Сравнение профилей температур вдоль канала

Среднеквадратичное отклонение профиля температур, полученного по ПК ЛОГОС и в эксперименте, составляет 0,053423. Рассчитанная температура в центрах твердотельных блоков отличается от результатов [35] не более чем на один градус, что соответствует относительной погрешности ~ 3%. Наблюдается достаточно хорошее качественное и количественное совпадение результатов.

Задача 20. Вентилируемая квадратная каверна с твердотельным блоком

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

Длина стороны каверны равна Ь , длина стороны твердотельного блока - ё. Внизу на левой стенке каверны расположен вход размером н'=0.1Ь, а наверху противоположной стенки находится выход такого же размера. Стенки каверны неподвижны. Верхняя, нижняя и левая вертикальная стенки теплоизолированные, а на правой вертикальной стенке поддерживается постоянная температура Тн. На входе в каверну задается температура Т, скорость входного потока щ направлена вдоль оси X . В данной задаче рассматривается вязкое ламинарное несжимаемое течение без учета вязкой диссипации. Учитывается вынужденная и естественная конвекция в приближении Буссинеска.

Рис. 68. Геометрия расчетной области

Постановка задачи взята из работы [36], где используется следующее обезразмерива-ние переменных:

X = £, Ь

у=У,

и = и,

г = ^,

р = -

о=±,

Ь

в =

т - т т - т

О, =

Т - т Ти - т

Ь иг и г Р"г

Здесь т и т - температура жидкости и твердотельного блока соответственно.

Число Рейнольдса Яе, число Ричардсона Яг, число Прандтля Рг и параметр К определяются следующим образом:

Яе =

иЬ

Яг =

- т)

Рг

к = $

где к„

у и2 к

и к - коэффициенты теплопроводности твердого тела и жидкости соответственно. Параметры жидкости выбирались таким образом, чтобы обеспечить значения чисел Яе=100, Рг=0.71 и К=5, указанные в работе [36]. Расчетная сетка регулярная ортогональная и равномерная, содержит 100 х100 х1 = 10000 ячеек.

В расчетах варьируется число Ричардсона Яг и размер твердотельного блока. Оценивается распределение температуры в жидкости и в твердом теле в сравнении с результатами работы [36], строится график зависимости числа Нуссельта от числа Ричардсона.

На рис. 69 представлены поля температуры, рассчитанные в ПК ЛОГОС при размере блока О = 0.4, с наложенными на них изолиниями температуры из работы [36]. На рис. 70 представлена аналогичная серия расчетов при О = 0.6. Видно практически полное совпадение изолиний, что свидетельствует о хорошем согласовании полей температуры.

Рис. 69. Рассчитанные поля температуры и изотермы при О = 0.4; Яг = 0, Яг = 1 и Яг = 5 соответственно Цветное поле - расчет ЛОГОС: черные изолинии -[36]

О ОООООООО О О О О О О О О ООО-'

О О^-*-—'■kjhowcjji bl Ch О? ОТ --j --j СО QO «D 4D О

О Ol О Ol О Ol О (Л О Ol О Ol О Ol О Ol О (Л О Ol о

Рис. 70. Рассчитанные поля температуры и изотермы из работы [36] при П = 0.6; Ш = 0,

Ш = 1 и Ш = 5 соответственно Цветное поле - расчет ЛОГОС, черные изолинии - [36]

Для численной оценки полученных результатов расчета проводится оценка числа Нуссельта на горячей стенке в зависимости от числа Ричардсона. Результаты представлены на рис. 71 и в табл. 16 и 17.

m 0 12.345

Ri

D = 0.4 D = 0.6

Рис. 71. Зависимость числа Нуссельта от числа Ричардсона

Таблица 16 Числа Нуссельта в зависимости от числа Ричардсона, П = 0.4

Ri 0,0 1,0 2,0 3,0 4,0 5,0

Nu, ЛОГОС 4,44 4,96 5,33 5,56 5,72 5,84

Nu, Rahman 4,3 4,85 5,24 5,47 5,63 5,75

S, % 3,3 2,3 1,7 1,6 1,6 1,6

Среднеквадратичное отклонение числа Нуссельта при размере блока П = 0.4 от результатов работы [36] составляет А = 0.0462.

Таблица 17 Числа Нуссельта в зависимости от числа Ричардсона, D = 0.6

Ri 0,0 1,0 2,0 3,0 4,0 5,0

Nu, ЛОГОС 4,71 4,94 5,17 5,36 5,54 5,68

Nu, Rahman 4,59 4,84 5,07 5,27 5,44 5,60

5, % 2,6 2,1 2,0 1,7 1,8 1,4

Среднеквадратичное отклонение числа Нуссельта при размере блока D = 0.6 от результатов работы [36] составляет А = 0.0443 .

Заключение

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

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

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

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

Библиографический список

1. AIAA, "Guide for the Verification and Validation of Computational Fluid Dynamics Simulations", AIAA G-077-1998, 1998.

2. Аксенова, Е.Н. Элементарные способы оценки погрешностей результатов прямых и косвенных измерений: учеб. пособие / Е.Н. Аксенова. - М.: Изд. МИФИ, 2003. - 16 с.

3. Ландау, Д. Гидродинамика / Д. Ландау, В.М. Лифшиц. - М.: Наука, 1988.

4. Козелков, А.С. Реализация метода расчета вязкой несжимаемой жидкости с использованием многосеточного метода на основе алгоритма SIMPLE в пакете программ ЛОГОС, ФГУП «РФЯЦ-ВНИИЭФ» / А.С. Козелков [и др.] // Журнал ВАНТ. Сер. Математическое моделирование физических процессов. 2013. Вып. 4. С. 44-56.

5. Голубев, А.А. Пакет программ ЛОГОС. Алгебраический многосеточный метод решения СЛАУ для задач гидродинамики / А.А. Голубев [и др.] // Современные проблемы науки и образования. 2013. № 6.

6. Швыдкий, В. С. Механика жидкости и газа / В. С. Швыдкий, Ю. Г. Ярошенко. - М.: ИКЦ «Академкнига», 2003. - 464 с.

7. Шлихтинг, Г. Теория пограничного слоя / Г. Шлихтинг. - М.: Наука, 1974.

8. Идельчик, И.Е. Справочник по гидравлическим сопротивлениям / И.Е. Идельчик. - М.: Машиностроение, 1992.

9. Козелков, А.С. Моделирование турбулентных течений вязкой несжимаемой жидкости на неструктурированных сетках с использованием модели отсоединенных вихрей / А.С. Козелков [и др.] // Матем. моделирование. 2014. 26:8 С. 81-96.

10.Ghia, U. High-Re Solution for Incompressible Flow Using the Navier-Stokes Equations and a Mul-tigrid Method / U. Ghia, K.N. Ghia, C.T. Shin // J. Comp. Phys. 1982. V. 48.

11.Stokes, G.G. Cambr. Trans. 1851. IX.

12.Mozer, D. &Mansour N. N. DNS of Turbulent Channel Flow / D. Mozer, J. Kim // Phys. Fluids. 1999. V. 11. Р. 943-945.

13.Ferziger, J.H. Computational Method for Fluid Dynamics / J.H. Ferziger, M. Peric // SpringerVerlag, New York, 2002.

14.Spalart, P.R. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach / P.R. Spalart [et al.] // Proceedings of first AFOSR international conference on DND/LES, 1997.

15.Spalart, P. R. Strategies for turbulence modeling and simulations // Heat Fluid Flow. 2000. V. 21. Р.252-263.

16.Travin, A. Physical and numerical upgrades in the detached-eddy simulation of complex turbulent flows / A. Travin [et al.] // Proceedings of Euromech Coll. Les of complex transitional and turbulent flows, Munich, Germany. Kluwer, Dordrecht. 2002. V. 65. Р. 239-254.

17.Spalart, P. R. A One-Equation Turbulence Model for Aerodynamic Flows. AIAA-92-0439.

18. Menter, F.R. Two-equation eddy viscosity turbulence models for engineering applications // AIAA J. 1994. 32. № 11. P. 1299-1310.

19.Menter, F.R. Ten Years of Industrial Experience with the SST Turbulence Model. Turbulence, Heat and Mass Transfer 4. Ed by K.Hanjalic, Y.Nagano, M.Tummers. Begell House Inc, 2003.

20.Бойко, А.В. Блок расчета положения ламинарно-турбулентного перехода для пакета программ ЛОГОС / А.В. Бойко [и др.] // Теплофизика и аэромеханика 2014. Т. 21. №2. С.201-220.

21.Лойцянский, Л.Г. Механика жидкости и газа / Л.Г. Лойцянский. - М., 1950.

22. База данных ERCOFTAG, тест 13 Source web site: ERCOFTAC Classic Database_Sudden Pipe Expansion (Experiments by Szczepura).

23.Iaccarino, G. Prediction of the turbulent flow in a diffuser with commercial CFD codes. Center for Turbulence Research, Annual Research Briefs 2000. Р. 271-278.

24.База данных ERCOFTAG, тест 80 Source web site: ERCOFTAC Classic Database_Steady Flow Past Tube Bundles (Experiments by S. Balabani).

25.Волков, К.Н. Моделирование крупных вихрей в расчетах турбулентных течений / К.Н. Волков, В.Н. Емельянов. - М.: Физматлит, 2008.

26. Jasak, H. Error Analysis and Estimation for the finite volume method with applications to fluid flows. Thesis submitted for the degree of doctor // Department of Mechanical Engineering, Imperial College of Science, 1996.

27. Leonard, B.P. A stable and accurate convective modeling procedure based on quadratic upstream interpolation // Comput. Methods Appl.Mech/Eng. 1979. V. 19. Р. 59-98.

28.Gaskell, P.H. Curvature-compensated convective-transport - SMART, A new boundedness-preserving transport algorithm // Int. J. Numer.Methods Fluids, 1988. V. 8. Р. 617-641.

29. Jasak, H. High resolution NVD differencing scheme for arbitrarily unstructured meshes / H. Jasak, H.G. Weller, A.D. Gosman // International journal for numerical methods in fluids, 1999. V. 31. Р.431-449.

30.Comte-Bellot, G. Simple Eulerian time correlation of full- and narrowband velocity signals in grid-generated "isotropic" turbulence / G. Comte-Bellot, S. Corrsin // Journal of Fluid Mechanics. 1971. V. 48. Р. 273-337.

31. Barakos, G. Natural convection flow in a square cavity revisited: laminar and turbulent models with wall fubctions / G. Barakos, E. Mitsoulis // Int. J. Numer. Methods Fluids. 1994. V. 18.

32. Wan, D.C. A new benchmark quality solution for the buoyancy-driven cavity by discrete singular convolution / D C. Wan, B.S. Patnaik, G.W. Wei // Numerical Heat Transfer, Part B, 40. 2001. Р.199-228.

33. Васильев, А.Ю. Инверсии крупномасштабной циркуляции при турбулентной конверсии в прямоугольных полостях / А.Ю.Васильев, П.Г. Фрик // Письма в ЖЭТФ. Т. 93. Вып. 6. С. 363-367.

34. Vogel, J.C. Combined Heart Transfer and Fluid Dynamic Measurements Downstream of a Backward-Facing Step / J.C. Vogel, J.K. Eaton // Journal of Heart Transfer. 1985. V. 107. С. 922-929.

35. Bilgen, E. Conjugate heat transfer in enclosures with openings for ventilation / E. Bilgen, T. Yama-ne // Heat and Mass Transfer 40. 2004. Р. 401-411.

36. Rahman, Md.M. Chowdhury, Effect of the presence of a heat conducting horizontal square block on mixed convection inside a vented square cavity / Md.M. Rahman [et al.] // Nonlinear Analysis: Modelling and Control. 2009. V. 14. № 4. Р 531-548.

37. Roache, P.J. Editorial Policy Statement on the Control of Numerical Accuracy / P.J. Roache [et al.] // ASME Journal of Fluids Engineering. Т. 108. 1986. № 1. С. 2.

Дата поступления в редакцию 08.10.2014

A.S. Kozelkov1'2, Yu.N. Deryugin1, Yu.A. Tsibereva1, A.V. Kornev3, O.V. Denisova1,

D.Yu. Strelets3, A.A. Kurkin 2, V.V. Kurulin1, I.L. Sharipova1, D.P. Rubtsova1, M.A. Legchanov2, E.S. Tyatyushkina1, S.V. Lashkin1, A.V. Yalozo1, S.V. Yatsevich1, N.V. Tarasova1, R.R. Giniyatullin1, M.A. Sizova1, O.L. Krutyakova1

MINIMAL BASIS TASKS FOR VALIDATION OF METHODS OF NUMERICAL SIMULATION OF TURBULENT FLOWS OF INCOMPRESSIBLE VISCOUS FLUIDS

FSUE «RFNC - VNIIEF»1, Nizhny Novgorod state technical university n.a. R.E. Alexeev , JSC "Sukhoi Company"3

Purpose: In this paper the validation process, which is an important step towards the industrial application of engineering software packages is discussed.

Method: The study is based on physico-mathematical and numerical models designed to simulate the turbulent flow of a viscous incompressible fluid.

Results: In order to improve and systematize the knowledge in this paper a basis of validation tasks required to assess the accuracy of the simulation software solutions designed for the simulation of turbulent flows of incompressible viscous fluid are developed. Held its systematization and generalization.

Application domain: Presented results allow the developer of engineering software packages focus on the computation of error modeling, rather than searching for reliable data.

Key words: validation, engineering software packages, turbulent flow, viscous incompressible fluid, numerical simulation.