Диссертация (1149258), страница 12
Текст из файла (страница 12)
Таким образом, интерес представляетоценка влияния методов ускорения RANS-LES перехода в слоях смешения на средниехарактеристики данного течения.Фотографии двухконтурного двигателя, использованного в экспериментах NASA,представлены на рисунке 3.10, где видны ротор, содержащий 22 лопатки, статор, онсостоит из 54 лопаток, а также центральное тело и кожух. Диаметр ротора составляет 22дюйма, а вся экспериментальная конструкция помещена в аэродинамическую трубуразмером 9x15 футов. Число Маха в трубе составляло M = 0.1, а скорость вращенияротора составляла 7808 об/мин.70Рисунок 3.10. Экспериментальная установкаВ рамках настоящей работы, как и в некоторых других работах ([107], [108],[110]), для расчетов использовалась упрощенная геометрия, включающая внутреннюючасть двигателя.
Помимо этого, для того чтобы уменьшить влияние входныхи выходных границ на течение, расчетная область включала в себя небольшую зонуперед кожухом и за ним (рис. 3.11).Рисунок 3.11. Вид полной геометрии двигателя. Красной линией отмечены контурыупрощенной геометрии, использованной в расчетах NASA [107], [108], а синей –использованной в настоящих расчетахРасчеты проводились в трехмерной постановке (геометрия и расчетная областьпредставлены на рисунке 3.12) в полной области (360 градусов).Вид расчетной сетки в плоскости XY, а также число ячеек в азимутальномнаправлении представлены на рисунке 3.12.
В области между ротором и статором шагсетки в осевом направлении не превышал 2.8∙10-3м, а общее число ячеек в направленииоси составило 480.71Рисунок 3.12. Трехмерный вид расчетной области, использованной в расчетах (сверху),а также расчетная сетка (снизу)На рисунке 3.13 представлен вид расчетной сетки в двух сечениях. Вблизи роторасетка сильно сгущена для того, чтобы разрешить течение в слоях смешения залопатками ротора, в то время как ниже по потоку распределение узлов в осевомнаправлении становится равномерным. В радиальном направлении во внутренней частидвигателя число ячеек расчетной области составляло 150.
При этом высота первойпристенной ячейки составляла 5∙10-6 м, что соответствует ее значению в переменныхзакона стенки менее 2 практически во всей области.Наконец, на рисунке 3.14 представлен вид расчетной сетки на поверхностилопаток статора и ротора, а также на центральном теле. Вокруг каждой лопаткистроился блок с сеткой О-типа.72Рисунок 3.13.
Сетка в плоскости YZ в двух сечения. Показана каждая вторая линиясеткиИспользованная в настоящих расчетах сетка близка к сетке, использованнойв работе [110] (упрощенная конфигурация) и состояла из 83 блоков. Общее число ячеексоставило порядка 55 миллионов.Рисунок 3.14. Расчетная сетка около лопаток ротора и статора, а также околоцентрального тела. Показана каждая вторая линия сеткиГраничные условия задавались следующим образом. На стенках задавалисьусловия прилипания и непроницаемости для скорости, а также условия адиабатичностидля температуры. На входной границе, расположенной примерно на расстоянии 0.3 мвыше по потоку от ротора задавался однородный профиль скорости. Радиальнаяи азимутальная компоненты вектора скорости полагались равными нулю.
Остальныепараметры течения на входной границе определялись с помощью характеристическихусловий с учетом заданного массового расхода Q = 26.44 кг/с и совпадения задаваемых73параметровсэкспериментальными(M∞ = 0.1,T∞ = 288 K,p∞ = 1.0135·105 Па).На выходной границе, расположенной на расстоянии 0.7 м от начала ротора задавалосьравномерное распределение статического давления, величина которого обеспечивалазаданный расход.
Остальные переменные на выходной границе определялисьс помощью линейной экстраполяции изнутри расчетной области.Вращение ротора обеспечивалось использованием вращающейся системыкоординат и заданием соответствующих граничных условий в части блоков (в блокахротораиобластимеждуроторомистатором).Награницевращающихсяи стационарных блоков использовалось условие интерфейса на основе одномернойинтерполяции (ширина перекрытия блоков составляла 3 ячейки).Шаг интегрирования по времени был равен 1/1056 части от времени полногооборота ротора.
Время выхода на установившийся режим составляло порядка 10 полныхоборотов ротора, а временная статистика собиралась за период, соответствующийвремени двух оборотов ротора.74Глава 4. Методы решенияКак уже отмечалось, результаты расчетов с помощью любых вихре-разрешающихподходов зависят от свойств используемых численных алгоритмов. Поэтому дляобъективнойоценкиэффективностиметодовускоренияRANS–LESперехода,получение которой является основной целью диссертации, для всех тестовых теченийкрайне желательно использовать единый численный метод.В настоящее время для аппроксимации невязких составляющих векторов потоковвисходныхуравненияхпереносаврамкахнезонныхгибридныхподходовобщепринятым является использование различных вариантов взвешенных центральноразностныхипротивопоточныхсхем,чтопозволяетодновременнодобитьсямонотонности решения в RANS подобласти, низкой диссипативности схемы в LESподобласти, и устойчивости алгоритма в целом.
Наиболее популярной из схем такоготипа является гибридная схема [17], в которой веса противопоточной и центральноразностной схем определяются автоматически на основе анализа текущего решенияи параметров используемой сетки. Однако эта схема была разработана для метода DES,применяемого, в первую очередь, для задач внешнего обтекания с массивнымиотрывными зонами, вследствие чего в пристенных течениях и течениях с умереннымиотрывными зонами она работает не оптимально: противопоточная схема зачастуюактивируется в пограничных слоях, населенных разрешенными турбулентнымиструктурами (LES подобласть), что приводит к диссипации разрешенных турбулентныхструктур, и в слоях смешения, что усугубляет проблему задержки RANS-LES перехода.В связи этим в рамках настоящей работы была разработана новая гибриднаячисленная схема для глобальных гибридных RANS-LES подходов, автоматическиобеспечивающая устойчивость вычислительного алгоритма в RANS подобласти и егонизкую диссипативность LES подобласти расчетной области и пригодная для болееширокого класса задач (разделы 4.1 и 4.2).
Эта схема была реализована на базевычислительного кода NTS [8] и использована для всех расчетов основной частиработы.NTS код – это конечнообъемный вычислительный код для решения уравненийНавье-Стокса, использующий структурированные многоблочные перекрывающиеся75расчетные сетки и подход, при котором основные переменные хранятся в узлах сетки(«vertex-basedapproach»).Кодпредназначендлярешениястационарныхи нестационарных течений как в сжимаемой, так и в несжимаемой постановкахс использованием различных моделей турбулентности. В рамках несжимаемой ветвикода применяется метод Роджерса и Квака [112], основанный на комбинации методавведения искусственной сжимаемости Яненко-Чорина [113] и схемы расщепленияразностей векторов газодинамических потоков.
Для расчетов сжимаемых теченийиспользуется схема Роу [114]. Для продвижения по времени используется неявнаятрехслойная схема Эйлера второго порядка точности. Для аппроксимации вязкихпотоков используется центрально-разностная схема второго порядка точности. Невязкиепотоки в рамках кода NTS могут аппроксимироваться с помощью противопоточной,центрально-разностной или взвешенной схемы разных порядков (от первого до пятого).Следует отметить, что код NTS рассматривается в настоящее время как один изнаиболее надежных вычислительных кодов для расчета турбулентных течений,полученные с его помощью результаты неоднократно сравнивались с результатамирасчетов, выполненных с использованием различных коммерческих и академическихкодов, предназначенных для решения задач гидро- и газодинамики (см., напримерработы [11], [12]).ВсерасчетыпроводилисьсиспользованиемгибриднойMPI-OpenMPпараллелизации на суперкомпьютере Политехнического университета.
Для различныхзадач использовалось от 1 до 7 узлов по 56 ядер каждый, а время расчета задачсоставляло от 1 до 20 дней.764.1. Формулировка новой гибридной численной схемы для глобальных гибридныхподходовВ рамках схемы [17] значения невязких потоков, аппроксимированныхс помощью центрально-разностной схемы (central-difference, CD) и противопоточнойсхемы, взвешиваются с помощью весовой функции, зависящей от расчетной сеткии мгновенного решения.Предложенная в настоящей работе новая схема также основана на взвесиневязких потоков, однако ее весовая функция отличается от использованной в работе[17].
Кроме того, в невязком потоке и в присоединенных пограничных слоях, в которыхотсутствует разрешенная турбулентность (RANS-область), эта схема переходитне в противопоточную схему, а во взвешенную схему BCD (Bounded CentralDifferencing, [115]), которая является менее диссипативной, чем противопоточная, но,тем не менее, обеспечивает монотонность решения и устойчивость алгоритма.
Принципработы схемы BCD основан использовании специального критерия («convectionboundedness criterion»), определяющего монотонность решения. В случае монотонностирешения используется центрально-разностная схема, в то время как в случаеневыполнения критерия, используется аппроксимация с помощью противопоточнойсхемы.Новая схема построена таким образом, чтобы во внешнем (невязком) потокеи в присоединенных пограничных слоях в отсутствие турбулентного контента дляаппроксимации невязких потоков использовалась схема BCD, в то время как в слояхсмешения, в том числе и на начальных их участках, а также в областях с разрешеннойтурбулентностью использовалась центрально-разностная схема либо взвесь центральноразностной схемы с небольшим весом BCD.Схема аппроксимации невязких потоков имеет вид:F (1 ) FCD FBCD .(4.1)Здесь FCD – значение невязких потоков, аппроксимированных с помощью центральноразностной схемы 4го порядка точности, FBCD – значение невязких потоков,аппроксимированныхспомощьюсхемыBCD,построеннойкаккомбинацияпротивопоточной схемы 3го порядка точности и центральной-разностной схемы 4гопорядка.
Наконец, σ – весовая функция, определяемая соотношением:77 max min , f inv , ft , f 2 D BL .(4.2)Здесь σmin – минимальный вес схемы BCD. Он принимает значения от 0 до 1 и задаетсяпользователем. В проведенных в настоящей работе расчетах использовалось значениепараметра σmin = 0. Функции finv, fνt, f2D BL используются для идентификации различныхобластей течения.Первая из перечисленных функций, finv, служит для определения невязкихобластей потока и определяется следующим образом: f inv 1 tanh B 4 ,(4.3)где B CH 3 max S , / max S 2 2 / 2, 2min- аргумент функции, предложеннойв работе [17], S и Ω – инварианты тензоров скоростей деформации и завихренностисоответственно, CH3 = 2.0, Ωmin = τ-1, а τ – характерное время конвективного проноса,построенное по характерным масштабам скорости и длины.