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

Полный текст

Введение. В настоящее время все большее значение придается исследованию проблем повышения надежности и ресурса жидкостных ракетных двигателей (ЖРД). В этой связи повышение надежности системы подачи компонентов топлива ЖРД, прежде всего турбонасосных агрегатов (ТНА), является ключевой проблемой. Узловым элементом ТНА является основной высокооборотный шнекоцентробежный насос компонента топлива (окислителя или горючего), обеспечивающий подачу рабочей жидкости в камеру сгорания и газогенератор при высоком давлении (свыше 100 бар). Шнекоцентробежный насос является основным источником гидродинамической вибрации системы подачи современных ЖРД. Гидродинамическая вибрация центробежного насоса является серьезной проблемой на пути повышения его надежности и ресурса на протяжении длительного времени. Первые упоминания об этой проблеме появляются в 60-х годах прошлого века в связи с разрушением крупных насосов [1; 2]. Гидродинамическая вибрация возбуждается пульсациями давления, возникающими в проточной полости насоса вследствие разных по своей природе гидродинамических причин [3; 4], которые включают вихреобразование, рециркуляцию потока, кавитацию, шаговую неравномерность параметров потока на выходе центробежного колеса. Последний фактор обусловливает генерацию пульсаций давления на так называемых лопаточных частотах или на частоте следования лопаток (ЧСЛ) и ее высших гармониках и комбинационных частотах. Эти колебания давления являются неотъемлемой частью рабочего процесса центробежного насоса. В центробежных насосах они имеют высокую амплитуду вследствие особенностей формирования шаговой неоднородности потока в центробежной лопаточной решетке. Хорошо известно, что физическая природа пульсаций давления в центробежном насосе представляет собой суммарное проявление псевдозвуковых и акустических колебаний. Псевдозвуковые колебания быстро затухают вниз по течению от ротора [5], оставляя в напорном трубопроводе только акустическую моду колебаний давления. Определение амплитуды пульсаций давления в шнекоцентробежном насосе на ранней стадии проектирования является актуальной задачей. В определении пульсаций давления, генерируемых трехмерным вихревым течением в шнекоцентробежном насосе, необходимо принимать во внимание, что неоднородное распределение параметров потока на выходе центробежного колеса генерирует акустические возмущения, которые распространяются со скоростью звука в рабочей жидкости. Одновременно присутствуют вихревые возмущения, которые конвектируются основным течением. Такие колебания параметров основного течения называют «псевдозвуком» [6; 7] или вихревой модой [8-10]. На двойственную природу пульсаций давления в центробежных насосах указывалось еще в ранних работах Покровского [3], далее указанный подход начал развиваться в работах Тимушева и соавторов [11] применительно к двумерной постановке задачи определения пульсаций давления в центробежном насосе с безлопаточным диффузором и спиральной улиткой. Течение в современных шнекоцентробежных насосах ЖРД имеет существенный трехмерный характер, поэтому возникает необходимость в адаптации акустико-вихревых уравнений для трехмерного случая. Адаптация акустико-вихревого метода для трехмерного течения. В данной работе метод расчета пульсаций давления с использованием подходов вычислительной гидродинамики и акустики развивается для случая трехмерного потока в шнекоцентробежном насосе с направляющим аппаратом и спиральным сборником. При этом используется представление нестационарного движения сплошной среды как совокупности акустической и вихревой мод движения. Физико-математическая модель генерации акустико-вихревых колебаний основана на представлениях, выдвинутых Ландау, Блохинцевым, Артамоновым [1-3], c применением декомпозиции граничных условий с использованием комплексного акустического импеданса [11]. Как уже отмечено выше, при pазpаботке физико-математической модели пульсационного течения в шнекоцентробежном насосе необходимо учитывать нелинейный характер процесса генерации колебаний нестационарным потоком и акустический характер их распространения в проточной части насоса. Приняты следующие допущения: - поток дозвуковой; - течение изоэнтропийное; - вязкая диффузия не учитывается для распространения акустических колебаний; - акустические колебания (вследствие сжимаемости среды) существенно меньше по сравнению с вихревыми колебаниями (вихревого и поступательного движения жидкости). Вводится акустическая и вихревая моды движения рабочей жидкости, так как на основании теоремы Коши-Гельмгольца скорость движения сжимаемой жидкости можно представить в виде векторной суммы основного поступательного и вращательного движения жидкости как несжимаемой среды (вихревой моды) и малых колебаний, обусловленных сжимаемостью среды (акустической моды). В качестве основной физической причины нестационарного процесса генерирования пульсаций давления в шнекоцентробежном насосе на частотах следования рабочих лопаток рассматривается перенос вихревых возмущений, которые возникают в результате движения периодически неоднородного потока с окружной скоростью U2 центробежного колеса относительно отводящего устройства. Такая модель генерации колебаний в литературе называется взаимодействием «ротор-статор». Затухание акустических возмущений, обусловленных вязкостью, а также тепловые явления имеют здесь второстепенное значение, и для упрощения соответствующие члены не будут учитываться в уравнениях движения, течение считается изоэнтропическим. Основные уравнения движения сжимаемой среды записываются в следующем виде: (1) (2) s = const. (3) (4) В изоэнтропическом течении приращения энтальпии, давления и плотности связаны термодинамическими соотношениями (5) где a - скорость звука в рабочей среде. С учетом соотношений (5), пренебрегая вязкостью в зоне распространения возмущений, перепишем уравнения (1) и (2) в форме, удобной для дальнейших преобразований: (6) (7) На основании теоремы Коши-Гельмгольца нестационарную скорость жидкости можно определить в виде суммы скорости U вихревой моды (поступательного и вращательного движения абсолютно несжимаемой среды) и скорости акустического движения Va. Введем скалярную функцию - акустический потенциал j. Тогда акустическая скорость (8) Таким образом, для скорости жидкости получается следующее выражение: (9) Метод разложения уравнений движения сплошной среды в форме (9) широко применяется в задачах аэроакустики [12] и гидродинамической устойчивости. В работе Артамонова [8], Тимушева [11] проведены преобразования основных уравнений движения в наиболее полном виде. Будем рассматривать дозвуковое течение M = = U/a << 1 с малыми акустическими колебаниями (Va << a). Запишем также следующие очевидные соотношения: (10) так как (11) Таким образом вихревые возмущения потока определяются скоростью несжимаемого течения. Подставим теперь соотношение (9) в уравнение (6). После несложных преобразований уравнение (6) приводится к виду (12) где (13) (14) В уравнении (12) член Ñj ´ x отражает взаимодействие акустической и вихревой мод. В этом уравнении учитывается действие вязких сил, которые участвуют в формировании неоднородности распределения параметров течения по шагу лопаточной решетки рабочего колеса. Теперь, выразив i из формулы (13) и подставив в уравнение (7), одновременно разделяя в нем скорость V на акустическую и вихревую моды, получим: (15) С учетом ранее принятых допущений можно линеаризовать уравнение (15) относительно j, и оно примет вид (16) Из уравнений (15) и (16) видно, что нестационарное вихревое движение жидкости генерирует акустические колебания. В качестве источника акустических колебаний в неоднородном волновом акустическом уравнении выступает в правой части уравнения член . В то же время энергия акустической моды может частично переходить в энергию вихревого движения, что выражается членом Ñj ´ x в уравнении (12). В данной модели член Ñj ´ x опускается, что существенно упрощает задачу численного решения трехмерных акустико-вихревых уравнений на компью-тере. С учетом последнего допущения и линеаризации по j соотношения (12), (13) записываются в следующем виде: (17) (18) Таким образом, уравнения (17), (18) обеспечивают решение задачи расщепления основных уравнений движения сжимаемой жидкости на вихревую и акустическую моды. Уравнение (17) описывает вихревое турбулентное движение несжимаемой вязкой жидкости под действием нестационарного градиента давления ÑPv = ρ0Ñj. Продифференцируем уравнение (15), взяв полную производную по времени, и подставим в него выражение для dj ¤ dt из формулы (18). После несложных преобразований получаем (19) Возмущающую функцию в правой части уравнения (19) можно выразить через поле скоростей вихревой моды из уравнения (17): (20) Пренебрегая конвективными членами [11] в производной по времени уравнения (19), получаем: (21) где через S обозначена возмущающая функция, определяемая из поля скоростей несжимаемого потока: (22) Перейдем к безразмерным переменным. В качестве пространственного масштаба и характерной скорости возьмем радиус R2 и окружную скорость u2 на выходе центробежного колеса. Тогда безразмерные величины запишутся в следующем виде: (23) Введем L - безразмерную частоту, или число Гельмгольца: (24) Как правило, в центробежных насосах L < 0,3, поэтому для низких гармоник амплитуд пульсаций давления ЧСЛ конвективные члены в волновом уравнении можно не учитывать: (25) Для невозмущенного потока акустический потенциал j = 0 и (26) Функции i и j можно выразить через средние величины и пульсационные составляющие: (27) Амплитуда пульсаций давления в отводе шнекоцентробежного насоса на порядок ниже среднего давления, поэтому для колебаний энтальпии можно приближенно записать (28) Аналогично для колебаний вихревой моды g (29) Учитывая соотношения (27)-(29), можно преобразовать формулу (18) к виду (30) Последнее выражение наглядно показывает, что пульсации давления рабочей жидкости в отводе шнекоцентробежного насоса равны сумме пульсаций от нестационарного вихревого движения (как несжимаемой среды) - «псевдозвука» и акустических колебаний. Принимая во внимание формулы (28), (29), из (25) получим уравнение первого приближения для колебаний давления: (31) где - нестационарная часть функции S. Для упрощения записи далее знак «%» в формулах опускается, и везде, где это не оговорено особо, рассматриваются безразмерные переменные. Решение уравнения (31) разделяется на две задачи: расчет нестационарного течения для модели несжимаемой среды, которое определяет возмущающую функцию, и решение неоднородного волнового уравнения относительно пульсаций давления h. Аналогичный подход используется, например, в работе [13]. Выражаем возмущающую функцию через поле скоростей несжимаемого течения из формулы: (32) (33) После преобразований на основе соотношений векторной алгебры получаем: (34) Учитывая уравнение неразрывности для несжимаемой среды, (35) Окончательно получаем следующее выражение: (36) Запишем акустико-вихревое уравнение в декартовой системе координат: (37) Возмущающая функция в правой части уравнения представляет собой нестационарную часть функции S (36). Конечно-разностные аналоги дифференциальных уравнений в декартовой системе координат получают-ся интегрированием акустико-волнового уравнения по пространству и времени с введением конечных объемов. Все поле течения покрывается прямоугольной сеткой. Каждому узлу сетки ставится в соответствие три числа (i, j, k), которые определяют порядковый номер конечного объема (ячейки) на Х-, Y- и Z-координатных осях. Границы между соседними ячейками проходят через середины шагов сетки. Кроме того, введем временную сетку с верхним индексом (m) и равномерным шагом по времени Dt, в которой каждому моменту времени соответствует номер m так, что (38) Рассмотрим вывод конечно-разностных уравнений для внутренних узлов сетки (ячейка 1). Для уравнения (37) используем интегральный метод [14]. Проинтегрируем уравнение (37) по пространству и времени в пределах одной ячейки и одного шага по времени: (39) здесь для внутренней ячейки 1 пределы интегрирования по объему (40) Будем считать, что давление, его вторая производная по времени и возмущающая функция постоянны внутри объема ячейки. Тогда можно записать для каждой гармоники ЧСЛ (41) С учетом формул для конечно-разностных аналогов производных (42) и уравнение (41) преобразуется к виду (43) где l - номер гармоники ЧСЛ; V - объем ячейки; функция Ψc (c - тип ячейки) определяется выражением (44) и зависит от типа ячейки, так как x1, y1, z1, x2, y2, z2 различны для разных типов ячеек. Для внутренних ячеек типа 1 (рис. 1) (45) Например, для ячеек, находящихся на стенке, нормальной к оси X справа, получаем (46) а для ячеек, находящихся на входной границе с нормалью n, (47) Таким образом, соответствующая производная от h по нормали к граничной поверхности определяется через производную по времени и удельный комплексный акустический импеданс Zi на данной границе для соответствующей гармоники ЧСЛ. Производные от колебаний вихревой моды известны из решения для нестационарного течения несжимаемой вязкой среды. Метод расчета вихревой моды. Расчет нестационарного течения в модели несжимаемой среды с использованием уравнения Навье-Стокса (48) и уравнения неразрывности (49) принят для первого шага акустико-вихревого метода: (48) (49) Эти уравнения дополнены уравнениями k-e (турбулентная энергия - скорость диссипации) модели турбулентности. В данной модели турбулентности [15] турбулентная вязкость mt выражается через величины k и e следующим образом: (50) Уравнения для k и e: (51) (52) где - начальное значение турбулентной диссипации; через G обозначено выражение (53) Значения параметров k-ε-модели равны: (54) в x2 - x1 = ∆x б x2 - x1 = ∆x а x2 - x1 = ∆x ris1т1@(2) Рис. 1. Тип ячейки: а - внутренняя; б - стенка; в - импедансная граница Граничное условие для скорости жидкости в турбулентном течении на стенке задается с использованием численной аппроксимации логарифмического закона для тангенциальной компоненты скорости на стенке [16]. Численный метод реализован на прямоугольной сетке с локальной адаптацией и подсеточным разрешением сложной геометрии. Во всей расчетной области вводится прямоугольная сетка. Выделяются подобласти с особенностями геометрии или течения, в которых необходимо провести расчет на более мелкой, чем исходная, сетке. При этом расчетная ячейка, в которую попала выделяемая особенность, делится на 8 равных ячеек. Далее, если необходимо, ячейки делятся еще раз, и так до достижения необходимой точности. Ячейки начальной сетки называются ячейками уровня 0, ячейки, получаемые измельчением уровня 0, называются ячейками уровня 1 и т. д. При генерации сетки накладывается условие, что гранями и ребрами могут граничить друг с другом только ячейки с номерами уровней, отличающимися не более чем на единицу. Метод подсеточного разрешения геометрии предназначен для аппроксимации криволинейных границ на прямоугольной сетке, в том числе и свободной границы жидкости. Ячейки, через которые проходит граница, расщепляются на 2, 3 и т. д. ячеек. При этом они теряют свою первоначальную форму параллелепипеда и превращаются в многогранники произвольной формы. Уравнения математической модели аппроксимируются для этих многогранников без каких-либо упрощений. Такой подход обеспечивает точность расчетов в зоне высоких градиентов параметров нестационарного потока (в зоне взаимодействия лопаточных решеток), а также эффективное решение волнового акустико-вихревого уравнения, используя минимальные вычислительные ресурсы. Геометрия проточной части и сетка. Для совмещения численных решений в роторе и статоре насоса применяется метод скользящих сеток. В этом случае вся область расчета разделяется на несколько подобластей, часть из которых решается во вращающейся системе отчета. Передача параметров потока из вращающейся в неподвижную область расчета производится через специальный интерфейс «скользящая поверхность», который обеспечивает интерполяцию параметров потока с учетом «виртуального» углового смещения сеток ротора и статора. В процессе решения акустико-вихревого волнового уравнения скользящие поверхности используются при постановке граничных импедансных условий для акустической моды пульсаций давления. Обмен данными между ротором и неподвижными подобластями выполняется с помощью интерполяции на границе подобластей - на «скользящих» поверхностях. Для этого временной шаг итерационной процедуры выбирается так, чтобы в течение одной итерации взаимное смещение подобластей не превышало размеров ячейки конечно-разностной сетки в зоне скользящих поверхностей. Расчетная область в математическом пространстве состоит из трех подобластей (рис. 2): подвода, ротора и статора, которые виртуально объединяются через скользящие поверхности, расположенные на входе и выходе шнекоцентробежного колеса и на ответных поверхностях подобластей подвода и статора. Данная расчетная область покрывается начальной прямоугольной сеткой, которая адаптируется к особенностям геометрии (рис. 3), обеспечивая необходимую точность расчета. Все три подобласти виртуально объединены по скользящим поверхностям. Таким образом, расчетная область адекватно отражает реальные условия течения рабочей жидкости, учитывая гидродинамическое взаимодействие всех элементов проточной части, включая направляющие ребра в подводе, вращающуюся лопаточную систему ротора с трехзаходным предвключенным шнеком и центробежным колесом, имеющим дополнительные укороченные лопатки (семь основных и семь дополнительных лопаток), а также двенадцатиканальный (трубчатый) направляющий аппарат. reg1 Рис. 2. Расчетная область grid_mer Рис. 3. Расчетная сетка второго уровня адаптации в меридиональной проекции Численное моделирование вихревой моды проводится нестационарным итерационным методом от нулевых начальных условий. Для выхода на стационарный (колебательный) режим течения в этом случае требуется обеспечить не менее шести полных оборотов центробежного колеса. Стационарный режим работы достигается для расчетной сетки при минимальной адаптации первого уровня с целью экономии процессорного времени и ресурсов. В качестве начальных условий в области жидкого объема задаются нулевые значения скорости и давления. Расчет ведется как процесс «раскрутки» насоса до достижения сходимости к периодическому осциллирующему решению. Как правило, шести-десяти полных оборотов рабочего колеса достаточно для получения периодического решения. Далее можно получить осреднением энергетические параметры насоса и сигналы и спектры пульсаций давления в разных точках проточной части направляющего аппарата. Сходимость численного решения контролируется по напору насоса. Временной шаг нестационарного расчета составляет 0,00001 с. Процессорное время расчета одного полного оборота ротора составляет около 2 ч. Утечки в уплотнениях колеса и сброс утечки на входе рабочего колеса не моделировались. Граничное условие заданного расхода применено в выходном сечении насоса. Последнее обеспечило возможность моделирования обратных токов на входе шнека при малых подачах. Статическое давление в расчетах отчитывается от опорной величины 101000 Па. Во входном сечении расчетной области задано статическое давление. Результаты численного моделирования. Проведено численное моделирование нестационарного течения в широком диапазоне изменения расхода относительно расчетного значения. При снижении расхода наблюдается перестройка течения в проточной части насоса как в подводе, так и в отводящем устройстве. В подводе на режимах 0,65 и особенно 0,3 относительного расхода фиксируются обратные токи на входе в шнек. На режиме 0,3 относительного расхода обратные токи проникают глубоко вверх по потоку, несмотря на наличие разделительного направляющего ребра. Такой характер течения вызывает дополнительные нестационарные нагрузки на корпус насоса и ротор. Кроме того, обратные токи вызывают круговую неоднородность потока в роторе, что приводит к усилению генерации колебаний на роторной частоте и на комбинационных частотах. Рассматривая основные тональные составляющие спектров в каналах направляющего аппарата, можно отметить, что существенное изменение амплитуды тональных компонент происходит при значениях относительного расхода 0,65-0,8, где при снижении расхода амплитуды основной лопаточной частоты и комбинационной гармоники на частоте 4 fr резко возрастают. При дальнейшем снижении расхода и усилении интенсивности обратных токов амплитуды лопаточных частот и комбинационной гармоники снижаются. В расчетах используется относительная амплитуда, приведенная к произведению плотности рабочей жидкости на квадрат окружной скорости на внешнем диаметре центробежного колеса. В каналах направляющего аппарата относительная амплитуда первой гармоники ЧСЛ (7 fr) на расчетном режиме составляет 0,0061. На рис. 4 показано изменение относительной амплитуды пульсаций давления первой гармоники ЧСЛ (7 fr) при изменении расхода по результатам акустико-вихревого моделирования на выходе насоса. На расчетном режиме расчетная величина относительной амплитуды составляет 0,00083. Это в 7 раз ниже, чем в направляющем аппарате. Известный экспериментальный уровень амплитуды пульсаций давления на выходе составляет 0,002. Такое несоответствие может быть связано с резонансным увеличением колебаний в реальном трубопроводе, так как в расчете использовано граничное условие акустического импеданса бесконечно длинной трубы. При относительном расходе от 0,8 до 1,35 от расчетного наблюдается плавный рост амплитуды. В зоне расхода менее 0,8 амплитуда пульсаций на данной частоте ведет себя нестабильно. Есть тренд к увеличению амплитуды при снижении расхода до 0,3 от расчетного. Такая нестабильность расчетных величин амплитуды, очевидно, связана с недостаточным осреднением источниковой функции акустико-вихревого уравнения: определение источниковой функции производится в течение одного периода движения лопаточной решетки. При наличии возмущений более низких частот, что характерно для режимов с обратными токами, точность определения амплитуды источниковой функции снижается. a1_vs_q Рис. 4. Изменение относительной амплитуды пульсаций давления на первой гармонике ЧСЛ [Па] в зависимости от относительного расхода через насос Заключение. Относительная амплитуда пульсаций давления ЧСЛ на выходе насоса (0,00083) в семь раз ниже, чем в каналах направляющего аппарата (0,0061). Значение амплитуды на выходе занижено по сравнению с известным экспериментальным уровнем амплитуды (0,002) вследствие влияния выходного акустического импеданса. При снижении относительного расхода до значений 0,65-0,8 амплитуда пульсаций давления ЧСЛ резко возрастает в направляющем аппарате и на выходе насоса.

Об авторах

С. Ф. Тимушев

Московский авиационный институт (Национальный исследовательский университет)

Российская Федерация, 125993, г. Москва, Волоколамское шоссе, 4

Д. В. Клименко

Московский авиационный институт (Национальный исследовательский университет)

Email: srmitriy@gmail.com
Российская Федерация, 125993, г. Москва, Волоколамское шоссе, 4

Список литературы

Дополнительные файлы