Решение выходного потока сверхзвукового сопла

Открыть скрипт

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

Описание задачи

В этом разделе описывается задача, которая будет решена. Это также предоставляет необходимые уравнения и известные значения.

Решение для поля потока после сверхзвукового сопла с помощью метода характеристик. Число Маха в плоскости выхода составляет 1,5, а давление в плоскости выхода составляет 200 килопаскалей. Заднее давление составляет 100 килопаскалей.

Допущения:

  • Поток изентропен

  • Изменение свойств потока зависит от взаимодействия волн расширения, которые происходят на протяжении следа сопла.

  • Геометрия сопла и потока симметричная.

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

upperNozzle = astexpandschematic('uppernozzle');

Приведенная информация в задаче:

exitMach = 1.5; % Mach number at the exit plane [dimensionless]
exitPres = 200; % Static pressure at the exit plane [kPa]
backPres = 100; % Pressure downstream of the nozzle, outside of the expansion wake

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

k = 1.4; % Specific heat ratio [dimensionless]

Метод характеристик

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

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

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

  • Угол потока является направлением, в котором воздух движется.

$$ \theta_n = Flow\;angle\;in\;the\;n^{th}\;region $$

  • Угол Прандтля-Мейера является углом, при котором поток изменяет направление от одной области к другой.

$$ \nu_n = Prandtl-Meyer\;angle\;in\;the\;n^{th}\;region $$

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

$$ \mu_n = Mach\;angle\;in\;the\;n^{th}\;region $$

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

Вычисление свойств потока через первый вентилятор расширения

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

[~, ~, exitPresRatio] = flowisentropic(k, exitMach);

Отношение противодавления является отношением противодавления к давлению застоя. Отношение изентропного давления во внешней области бодрствования составляет:

$$ \frac{p_4}{p_0} = \frac{p_4}{p_1} \frac{p_1}{p_0}$$

backPresRatio = backPres / exitPres * exitPresRatio;

Вычислим число Маха в области 4 с помощью flowisentropic.

backMach = flowisentropic(k, backPresRatio, 'pres');

Вход строки 'pres' указывает, что функция находится в режиме входа отношения давления. Угол потока в области заднего давления является различием в углах Прандтля-Мейера от области выходной плоскости (область 1) до области заднего давления (область 4).

$$ \theta_4 = \nu_4 - \nu_1 $$

[~, nu_1]    = flowprandtlmeyer(k, exitMach);
[~, nu_4]    = flowprandtlmeyer(k, backMach);
theta_4 = nu_4 - nu_1;

Поскольку мы аппроксимируем поток с тремя характеристиками, вычислите изменение угла потока при пересечении характеристик типа I от области 1 до области 4:

$$ \Delta \theta_I = \frac{\theta_4}{3} $$

deltaThetaI = theta_4 / 3;

Обратите внимание, что поток в области 1 параллелен горизонтали и поэтому:

theta_1 = 0;

Фактически, поток в любой области, которая соединяет осевую линию, параллелен осевой линии. Это связано с тем, что осевая линия рассматривается как контур для этого симметричного потока. В сложение на контур нет ни источников, ни стоков.

theta_5  = 0;
theta_8  = 0;
theta_10 = 0;

Углы потока областей 2 и 3 следуют просто.

theta_2 = theta_1 + deltaThetaI;
theta_3 = theta_2 + deltaThetaI;

Между характеристиками типа I изменение угла Прандтля-Мейера равно изменению угла потока:

$$ \Delta \nu_I = \Delta \theta_I $$

deltaNuI = deltaThetaI;

Вычислите угол Прандтля-Мейера в области 2 с помощью угла Прандтля-Мейера в области 1 и deltaNuI, изменения угла Прандтля-Мейера через характеристики типа I. Вычислите угол Прандтля-Мейера в области 3 так же, как в области 2.

$$ \Delta \nu_I = \nu_2 - \nu_1 $$

$$ \nu_2 = \nu_1 + \Delta \nu_I $$

nu_2 = nu_1 + deltaNuI;
nu_3 = nu_2 + deltaNuI;

Вычисление свойств потока в интерференционных областях

Известно, что угол потока в области 5 равен нулю от граничного условия осевой линии. Поэтому изменение угла от области 2 к области 5

$$ \Delta \theta_{II} = \theta_5 - \theta_2 $$

deltaThetaII = theta_5 - theta_2;

Вычислите изменение угла Прандтля-Мейера между характеристиками типа II:

$$ \Delta \nu_{II} = - \Delta \theta_{II} $$

deltaNuII = -deltaThetaII;

Затем вычислите угол Прандтля-Мейера в области 5. Вы уже знаете область 2 углом Прандтля-Мейера и изменением угла Прандтля-Мейера между характеристиками типа II.

nu_5 = nu_2 + deltaNuII;

Для вычисления свойств в области 6 используйте тот факт, что свойства в области 3 и области 5 известны. Обратите внимание, что характеристики пересечения потока определяют изменения свойств. От области 5 до области 6 пересекается характеристика I типа. Поэтому,

$$ \Delta \theta_I = \Delta \nu_I = \theta_6 - \theta_5 = \nu_6 - \nu_5 $$

Переставили это как:

$$ \nu_6 - \theta_6 = \nu_5 - \theta_5 \;\;\;\;\;(1)$$

Характеристику типа II пересекают при переходе от области 3 к области 6. Поэтому,

$$ \Delta \nu_{II} = -\Delta \theta_{II} $$

$$ \nu_6 - \nu_3 = \theta_3 - \theta_6 $$

Реорганизация этого параметра как:

$$ \nu_6 + \theta_6 = \nu_3 + \theta_3 \;\;\;\;\;(2)$$

Добавьте уравнения (1) и (2) вместе, затем решите для угла Прандтля-Мейера в области 6. Это дает следующее выражение.

$$ \nu_6 = \frac{(\nu_5 - \theta_5) + (\nu_3 + \theta_3)}{2}$$

В MATLAB ® используйте:

nu_6 = ((nu_5 - theta_5) + (nu_3 + theta_3))/2;

Из уравнения (1) угол потока в области 6 равен

theta_6 = nu_6 - (nu_5 - theta_5);

Для области 7 пересекается характеристика типа один, и вся информация доступна в области 6.

nu_7    = nu_6 + deltaNuI;
theta_7 = theta_6 + deltaThetaI;

Область 8 расположена на осевой линии; его угол потока равен нулю. Переход от области 6 к области 8 требует пересечения характеристики типа II. Поэтому вычислите угол Прандтля-Мейера в области 8 как:

nu_8 = nu_6 + deltaNuII;

Вычислите угол Прандтля-Мейера и угол потока в области 9 так же, как вы это сделали для области 6. Область 8 является верхней областью через характеристику типа I. Область 7 является верхней областью через характеристику типа II.

nu_9    = ((nu_8 - theta_8) + (nu_7 + theta_7))/2;
theta_9 = nu_9 - (nu_8 - theta_8);

Область 10 расположена на осевой линии. Поток параллелен, и поэтому угол потока равен нулю. Используйте угол Прандтля-Мейера в области 9 и пересечение по характеристике типа II, чтобы вычислить угол Прандтля-Мейера в области 10.

nu_10 = nu_9 + deltaNuII;

Подготовка и табуляция результатов параметра потока

Для предстоящих вычислений объедините углы потока в один вектор и углы Прандтля-Мейера в другой вектор.

flowAngles         = [theta_1 theta_2 theta_3 theta_4 theta_5 theta_6 theta_7 theta_8 theta_9 theta_10];

prandtlMeyerAngles = [nu_1 nu_2 nu_3 nu_4 nu_5 nu_6 nu_7 nu_8 nu_9 nu_10];

Чтобы вычислить числа Маха и углы Маха в каждой области, используя flowprandtlmeyer функция с prandtlMeyerAngle в качестве входов. Можно использовать результаты этой функции, чтобы найти угол, который имеют признаки типа I и типа II с горизонтальными внутри каждой области. Затем можно использовать эти углы для вычисления уклонов в плоскости x-y, где осевая линия является осью X, а выходная плоскость сопла - осью Y. Для характеристик типа I и типа II, соответственно, уклоны:

$$ \left(\!\frac{dy}{dx}\!\right)_{\! I} = tan(\theta - \mu)$$

$$ \left(\!\frac{dy}{dx}\!\right)_{\!\! II} = tan(\theta + \mu)$$

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

% Preallocation for speed
machNumbers = zeros(1,10);
machAngles  = zeros(1,10);
typeOne     = zeros(1,10);
typeTwo     = zeros(1,10);

for i = 1:10
    [machNumbers(i), ~, machAngles(i)] = flowprandtlmeyer(k, prandtlMeyerAngles(i), 'nu');
    typeOne(i) = flowAngles(i) - machAngles(i);
    typeTwo(i) = flowAngles(i) + machAngles(i);
end


clear table;
table(1,:) = 'Region   theta        nu        Mach        mu       type I    type II';
table(2,:) = '         (Deg)      (Deg)                  (Deg)     (Deg)      (Deg) ';
for m=1:length(machNumbers)
table(m+3,:) = sprintf('%3.0d   %8.2f      %5.2f   %8.3f    %8.2f   %8.2f  %8.2f ', ...
                m, flowAngles(m), prandtlMeyerAngles(m), machNumbers(m), ...
                machAngles(m), typeOne(m), typeTwo(m));
end

disp(table)

Region   theta        nu        Mach        mu       type I    type II
         (Deg)      (Deg)                  (Deg)     (Deg)      (Deg) 
                                                                      
  1       0.00      11.91      1.500       41.81     -41.81     41.81 
  2       4.45      16.35      1.650       37.29     -32.85     41.74 
  3       8.89      20.80      1.803       33.70     -24.80     42.59 
  4      13.34      25.24      1.959       30.69     -17.35     44.03 
  5       0.00      20.80      1.803       33.70     -33.70     33.70 
  6       4.45      25.24      1.959       30.69     -26.25     35.14 
  7       8.89      29.69      2.122       28.11     -19.22     37.00 
  8       0.00      29.69      2.122       28.11     -28.11     28.11 
  9       4.45      34.14      2.294       25.84     -21.40     30.29 
 10       0.00      38.58      2.477       23.81     -23.81     23.81

Обратите внимание на следующее:

  • Углы потока увеличиваются от осевой линии.

  • Углы Прандтля-Мейера увеличиваются, когда поток движется вниз по течению.

  • Число Маха также увеличивается, когда поток движется вниз по потоку.

Решение для геометрии потока

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

Предположим, что пересечение осевой линии и выходной плоскости сопла является источником нашей системы координат. Также предположим, что длины нормированы к половине высоты выхода сопла. Положительная ось x принята горизонтальной вдоль осевой линии в направлении вниз по течению. Положительная ось y расположена вертикально вверх в выходной плоскости сопла. Выступ сопла находится в точке (0,1).

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

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

Самая крутая линия наклона от выступа является характеристикой типа I, которая разделяет область 1 и область 2. Чтобы вычислить угол, который самая крутая наклонная волна делает с горизонтальной используйте среднее значение углов, которые характеристики типа I в каждой области. Для вычисления уклона используйте тригонометрию.

avgAngle12 = (typeOne(1) + typeOne(2)) / 2;
slope12    = tand(avgAngle12);

Со следующей известной информацией:

  • Наклон первого типа I волны в пространстве x-y.

  • Свободный член волны (y = 1 на губе).

  • Волна пересекает осевую линию (y = 0) без пересечений.

Вычислите положение точки по оси X с помощью уравнения линии в форме «уклон-перехват». Переставьте y = m * x + b для местоположения x с y = 0, чтобы получить x = -b/m. Это местоположение x первой нисходящей точки, точка 1.

y1 = 0; % On the centerline
x1 = -1 / slope12;

С точки 1 первая характеристика типа II распространяется и мешает вентилятору. Другие характеристики типа I, которые берут начало от выступа сопла, нарушаются волной типа II, но не до достижения этой волны. Поэтому вычислите точки пересечения самой крутой характеристики типа II и более плоских волн типа I от губы. Характеристика типа II, выходящая из осевой линии, разделяет область 2 и область 5. Среднее значение двух углов и связанных склонов определяется:

avgAngle25 = (typeTwo(2) + typeTwo(5)) / 2;
slope25    = tand(avgAngle25);

Вторая самая крутая характеристика типа I является от выступа сопла, разделяет область 2 и область 3. Средний угол с горизонтальным и связанным с ним наклоном волны задается:

avgAngle23 = (typeOne(2) + typeOne(3)) / 2;
slope23    = tand(avgAngle23);

Вычислите точку пересечения области 2-3 контуров и области 2-5 контуров. Вам нужна эта точка, потому что характеристики мешают друг другу на данном этапе. Известны склоны обоих контуров и точки на каждой линии. Точка 1 и выступ сопла (для ссылки на точку 0 известны ". Решение для неизвестной координаты Х точки пересечения. Используйте это положение X в уравнении любой из двух линий, чтобы найти положение y точки пересечения. Форма точки-уклона уравнения линии через точку p с уклоном m является:

$$ y - y_p = m(x - x_p)$$

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

$$ y_2 - y_0 = m_{2,3}(x_2 - x_0)\;\;\;\;(3)$$

$$ y_2 - y_1 = m_{2,5}(x_2 - x_1)\;\;\;\;(4)$$

где

$$ m_{i,j} = Slope\;of\;the\;boundary\;between\;region\;i\;and\;region\;j$$

Вычесть и переставить:

$$ x_2 = \frac{x_1 m_{2,5} - x_0 m_{2,3} + y_0 - y_1}{m_{2,5} - m_{2,3}}$$

Знание некоторых значений точно из-за точек пересечения осей упрощает это выражение в.

x2 = (x1 * slope25 + 1) / (slope25 -slope23);

Ниже y-местоположения точки 2 определяется путем закупоривания x-расположения точки 2 в уравнение (4) выше, но закупоривание в уравнение (3) работает также.

y2 = (x2 - x1) * slope25;

Используйте формулу перехвата наклона и процедуру выше, чтобы вычислить все точки. Чтобы вычислить третью точку потока, сначала вычислите пересечение области 3-4 контуры и 3-6 контуры. Углы граничных линий вычисляются с помощью среднего значения углов с горизонтальными. Затем можно использовать тригонометрию, чтобы найти уклон, который теперь вычисляется за один шаг.

slope34 = tand( (typeOne(3) + typeOne(4)) / 2 );
slope36 = tand( (typeTwo(3) + typeTwo(6)) / 2 );

Поскольку контур между областью 3 и областью 4 является характеристикой типа I, а контур между областью 3 и областью 6 является характеристикой типа II, будьте осторожны, чтобы взять углы для соответствующего типа. Используйте форму «точка-уклон» этих граничных линий, вычитаемых друг из друга, чтобы вычислить положение точки 3 по оси X.

x3 = (y2 - 1 - x2 * slope36) / (slope34 - slope36);
y3 = (x3 - x2) * slope36 + y2;

Первая характеристика типа II теперь распространяется за пределы вентилятора и не препятствует каким-либо другим характеристикам. Угол, который определяет направление, в котором распространяется самый крутой тип II, является углом, который является средним значением волн типа II в областях 4 и 7.

slope47 = tand( (typeTwo(4) + typeTwo(7)) / 2);

Решите также вторую самую крутую характеристику типа I, которая будет продолжаться ниже по потоку. Начните с известного местоположения точки 2, чтобы вычислить x-перехват второй характеристики типа I. Снова в решении используется форма «точка-уклон» линии. Однако пересечение отсутствует до достижения контура осевой линии. Нам нужно рассмотреть только одну линию, чтобы найти местоположение X «точки 4». Наклон контура между областью 5 и областью 6 является волной I типа:

slope56 = tand( (typeOne(5) + typeOne(6)) / 2 );

Переставленная форма уклон-точка (зная y = 0 на осевой линии) используется для поиска точки 4.

x4 = ( slope56 * x2 - y2 ) / slope56;
y4 = 0;

Вычислите точку 5 так же, как и точку 2. Область 6-7 контуров является типом I, а область 6-8 контуров - типом II.

slope67 = tand( ( typeOne(6) + typeOne(7)) / 2);
slope68 = tand( (typeTwo(6) + typeTwo(8)) / 2);

Известная точка на контуре области 6-7 является точкой 3. Известная точка на контуре области 6-8 является точкой 4. Используйте эту информацию в форме slope-intercept, вычитая уравнения, и переставьте, чтобы получить местоположение следующей точки.

x5 = (-x4 * slope68 + x3 * slope67 + y4 -y3) / (slope67 - slope68);
y5 = (x5 - x4) * slope68 + y4;

Вторая характеристика типа II распространяется за пределы вентилятора под углом, усредненным между областью 7 и областью 9 углами для волн типа II.

slope79 = tand( (typeTwo(7) + typeTwo(9)) / 2);

Последней точкой интереса является x-перехват самой плоской волны I типа. Вычислите эту точку путем знания местоположения точки 5 и нахождения наклона волны типа I между областью 8 и областью 9.

slope89 = tand( (typeOne(8) + typeOne(9)) / 2);
y6      = 0;
x6      = (slope89 * x5 - y5) / slope89;

Конечная волна II типа распространяется под углом, усредненным между областью 9 и областью 10.

slope910 = tand( (typeTwo(9) + typeTwo(10)) /2);

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

points = astexpandschematic('nozzlepoints');

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

close(upperNozzle,points)

Ссылка

[1] James, J. E. A., «Gas Dynamics, Second Edition», Allyn and Bacon, Inc, Boston, 1984.

Документация Aerospace Toolbox

Поддержка

Для просмотра документации необходимо авторизоваться на сайте