Журнал вычислительной математики и математической физики, 2021, T. 61, № 10, стр. 1684-1692

ВВЕДЕНИЕ

С каждым годом возрастает уровень требований к характеристикам летательных аппаратов. Поэтому задачи оптимизации элементов летательных аппаратов и их двигателей, позволяющие улучшать эти характеристики, еще долгое время будут оставаться актуальными. В работе рассматривается задача о нахождении формы пространственной сверхзвуковой части сопла, вписанной в заданные габариты и имеющей наибольшую тягу среди всех возможных допустимых форм. В настоящее время применяются два подхода к решению этой задачи: прямые приближенные методы и обратные точные методы. В прямых методах искомая стенка параметризуется, и искомые параметры находятся с помощью прямых методов максимизации интеграла тяги (см. [1]–[3]). Обратные методы используют математический аппарат теории оптимального управления системами с распределенными параметрами, когда связи задачи учитываются с помощью множителей Лагранжа. Далее вычисляется первая вариация расширенного функционала, возникают сопряженная система для множителей Лагранжа и условие оптимальности. Совместное решение исходной и сопряженной систем, удовлетворяющее условию оптимальности, позволяет найти оптимальное решение. Но решение сопряженной системы для множителей Лагранжа с достаточной точностью является весьма трудоемким процессом, который стал возможным только теперь с развитием вычислительной техники. В этой процедуре для задач с распределенными параметрами есть “узкое” место: какие связи задачи необходимо учитывать с помощью множителей Лагранжа, а какие можно и не привлекать? Дело в том, что учет всех связей в применимости к пространственным оптимальным задачам приводит к чудовищно громоздким расширенным функционалам и громоздким сопряженным системам. Авторы монографий Ю.Д. Шмыглевский (см. [4]) и А.Н. Крайко (см. [5]) считают, что необходимо учитывать только существенные связи задачи. В некоторых задачах можно получить оптимальные решения, которые автоматически удовлетворяют не учтенным в функционале Лагранжа связям. Последнее связано с тем, что если бы такие связи были учтены с помощью соответствующих множителей Лагранжа, то эти множители все равно получились бы нулями (см. [4]). Примером такой задачи является задача, рассматриваемая в настоящей статье. В ней не учитываются условия совместности на характеристических поверхностях, ограничивающих область решения.

1. ПОСТАНОВКА ЗАДАЧИ

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

Система (1) для идеального газа замыкается интегралом Бернулли и уравнением состояния где $\unicode{230} $ – показатель адиабаты.

Все искомые функции, входящие в уравнения (1)–(3): скорость V, давление p, плотность ρ, а также местная скорость звука $a = \sqrt {\unicode{230} \frac{p}{\rho }} $ обезразмерены к параметрам в критическом сечении сопла. Независимые переменные отнесены к радиусу критического сечения.

Из формул (1)–(3) следует, что закон сохранения импульса в проекции на фиксированное направление ${\mathbf{e}}$ можно записать в следующем дивергентном виде:

Дополнительно введем две скалярные функции тока ψ и χ, определяемые уравнениями

Функции ψ и χ сохраняются вдоль линий тока, поэтому задавая значения функций тока ψ и χ на некоторой начальной поверхности, пересекающей все линии тока, мы можем с помощью (5) рассчитать значения ψ и χ во всей области течения и выделить нужные линии тока.

Начальные данные задаются на некоторой характеристической поверхности Σ₁, сходящей с круглого критического сечения сопла γ₁ (фиг. 1). На стенке сверхзвуковой части сопла σ выполняется условие непротекания

где n_σ – нормаль к поверхности σ.

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

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

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

Расчет функций V, p, ρ, ψ и χ в заданной области τ (фиг. 1), ограниченной начальной характеристической поверхностью Σ₁, замыкающей характеристической поверхностью Σ₂, приходящей на заданный выходной контур сопла γ₂, и заданной стенкой сверхзвуковой части сопла σ, по уравнениям (1)–(3), (5) с начальными данными на характеристической поверхности Σ₁ будем называть прямой задачей.

Сформулируем теперь вариационную задачу для определения оптимальной формы сверхзвуковой части сопла. Пусть заданы входная критическая γ₁ и выходная γ₂ кромки сверхзвуковой части сопла (в общем случае контур γ₂ может лежать на известной поверхности Ω, см. фиг. 1). Пусть далее вектор e задает направление, в котором будет проводиться оптимизация тяги сопла. Требуется среди всех поверхностей σ, проходящих через заданные контуры γ₁ и γ₂, найти такую, на которой функционал

являющийся проекцией интеграла сил давления на направление e, принимает наибольшее значение (здесь p₀ – постоянное атмосферное давление на внешней поверхности сопла). Давление p в (8) находится из решения задачи (1)–(3). В точках критического сечения продольные образующие сопла имеют излом.

2. СОПРЯЖЕННАЯ ЗАДАЧА

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

в котором дифференциальные уравнения (1) и условие непротекания (6) учитываются с помощью переменных множителей Лагранжа $~{{\mu }}$, λ и h. Тогда задача нахождения экстремума функционала (8) при условии выполнения уравнений (1)–(3) сводится к задаче нахождения безусловного экстремума функционала T⁰. Таким образом, необходимым условием экстремальности сверхзвуковой части сопла является обращение первой вариации функционала δT⁰ в нуль.

Для получения первой вариации мы будем использовать следующие формулы, основанные на обобщении результатов Н.Е. Кочина по переменным полям в сплошных средах (см. [7]). Формулы, совпадающие с формулами Н.Е. Кочина, опубликованы также и в монографии [8]. Пусть заданы функционалы

Тогда первые вариации этих функционалов равны Здесь вектор-функция δa обозначает допустимую вариацию вектора a внутри и на границе фиксированной поверхности S, δφ – допустимая вариация функции φ внутри и на границах фиксированного объема τ, вектор-функция δr_S – малая непрерывная деформация поверхности S, ограничивающей объем τ, δr_L – малая непрерывная деформация контура L, ограничивающего поверхность S.

При варьировании поверхности σ из-за сверхзвукового характера течения будем считать, что на поверхности Σ₁ функции течения не меняются. Это означает, что на поверхности Σ₁ выполняются соотношения $\delta {\mathbf{V}} = \delta {\mathbf{r}} = 0$. Кроме того, из формул (1)–(3) следуют связи

С учетом этих формул первая вариация функционала Лагранжа записывается в виде где n₂ – вектор нормали к поверхности Σ₂ такой, что ${\mathbf{V}} \cdot {{{\mathbf{n}}}_{2}}~$ = a.

Определим множители Лагранжа следующим образом:

в области τ

на поверхности Σ₂

на поверхности σ

Выпишем теперь три дифференциальных и одно конечное следствия формул (9)–(11), не содержащих вектор-функцию h:

в области τ

на поверхности Σ₂

на поверхности σ

на контуре γ₂

Уравнения (13)–(16) однозначно определяют функцию λ в замкнутой области τ. Дело в том, что уравнение (13) имеет гиперболический тип и его характеристические поверхности совпадают с характеристическими поверхностями уравнений (1). Значение нормальной производной ${{{\mathbf{n}}}_{\sigma }} \cdot \nabla \lambda $ известно на σ$~$из (15), сама функция $\lambda ({{{{\Sigma }}}_{2}})$ однозначно определяется уравнением (13) и граничными условиями (14), (16).

Отметим, что после определения множителя $\lambda (\sigma )$ формула (12) однозначно определяет множитель $\mu (\sigma )$.

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

Теперь условие оптимальности поверхности σ можно записать с учетом формулы (4) в виде Интегралом последнего равенства является соотношение где $\Phi (\psi ,\chi )$ – произвольная функция, сохраняющая свое значение на линиях тока.

Вектор dr во втором интеграле (17) лежит на поверхности Ω, поэтому его можно представить в точках контура γ₂ в виде

где N – вектор нормали к поверхности Ω.

Окончательно получим

Если контур γ₂ фиксирован, то на нем $\delta {{{\mathbf{r}}}_{L}} = 0$, если контур γ₂ не фиксирован и принадлежит поверхности Ω, то вдоль кромки γ₂ с учетом (16) получаем аналог локального условия Буземана для пространственных вариационных задач газовой динамики

В вариационных задачах о соплах максимальной тяги часто существуют конструктивные ограничения на контур γ₂, приводящие к “недорасширенным” соплам. В таких соплах габаритные ограничения на кромку γ₂ не позволяют получить степень расширения, соответствующую условию (20). Поэтому на участках краевого экстремума допустимая в формуле (19) произвольная вдоль γ₂ функция ${{{\mathbf{n}}}_{\sigma }} \cdot \delta {{{\mathbf{r}}}_{L}} < 0$. Отсюда на таких участках необходимое условие локального максимума тяги вместо равенства (20) приводит к неравенству

В дальнейшем приведем примеры расчетов только оптимальных “недорасширенных” сопел, вдоль выходного контура которых выполняется достаточное условие локального максимума: строгое неравенство в формуле (21).

Введем новую искомую функцию ${{\lambda }_{1}} = \lambda + {\mathbf{e}} \cdot {\mathbf{V}}$. Тогда формулы (13)–(16) можно переписать в виде

в области τ

на поверхности Σ₂

на поверхности σ

на контуре γ₂

Условие оптимальности (18) на поверхности сопла $\sigma $ примет вид

В дальнейшем систему (22)–(25) будем называть сопряженной задачей.

Из уравнения (22) видно, что вектор $\rho \nabla {{\lambda }_{1}} - \frac{\rho }{{{{a}^{2}}}}{\mathbf{V}}({\mathbf{V}} \cdot \nabla {{\lambda }_{1}})$ является соленоидальным и, следовательно, может быть представлен как вихрь некоторого другого вектора A:

Вектор A в общем случае можно представить в виде трех скалярных функций, т.е. мы имеем всего три скалярных уравнения (27) для определения четырех неизвестных функций λ₁, A. Поэтому однозначно определить эти функции из (27) мы не можем. Воспользуемся этим произволом и представим вектор A с помощью двух неизвестных функций λ₂ и λ₃. Такое представление неединственно. Возьмем его в виде Здесь φ – угол цилиндрической системы координат x, r, φ (орты ${{{\mathbf{e}}}_{x}},\;{{{\mathbf{e}}}_{r}},\;{{{\mathbf{e}}}_{\varphi }}$), в которой мы и будем работать в дальнейшем. Выбор такого представления вектора A можно объяснить двумя причинами. Во-первых, в осесимметричном случае при λ₃ = 0 уравнения (27) совпадают с уравнениями, полученными в [6]. Во-вторых, как это будет видно дальше, краевые условия для λ₂ и λ₃ имеют наиболее простой вид. Таким образом, система уравнений (27) записывается в виде

Система уравнений (28) имеет гиперболический тип. Ее характеристические поверхности совпадают с характеристическими поверхностями уравнений (1). Условие совместности на них получается скалярным умножением (28) на вектор нормали к характеристической поверхности. В частности, условие совместности на замыкающей характеристической поверхности Σ₂ имеет вид

Выберем распределения функций ψ и χ в критическом сечении сопла следующим образом: Тогда в точках контура γ₁ и, следовательно, на всей стенке сопла ψ = 1. Умножим скалярно уравнение (28) на вектор n_σ. Так как $\sigma $ является поверхностью уровня для функции ψ, то $\nabla \psi \times {{{\mathbf{n}}}_{\sigma }}$ = 0 и $~{\mathbf{V}} \cdot {{{\mathbf{n}}}_{\sigma }}$ = 0, и второе слагаемое в правой части (28) пропадает. Слагаемое в левой части в силу (24) также обращается в нуль. Окончательно получим условие, справедливое на стенке сопла, $({{{\mathbf{e}}}_{\varphi }} \times {{{\mathbf{n}}}_{\sigma }}) \cdot \nabla {{\lambda }_{2}} = 0$. Отсюда следует, что на поверхности сопла мы имеем условие λ₂ = const вдоль линий пересечения σ с меридиональными плоскостями. И если мы зададим λ₂ = 0 во всех точках контура γ₂, то на всей поверхности σ Если теперь задать на всей поверхности Σ₂ какое-либо распределение λ₃, например, λ₃ = 0, то из уравнений (23), (29) с граничными условиями (24), (25), (30) можно найти распределения λ₁ и λ₂ на поверхности Σ₂. Таким образом, в качестве краевых условий для системы (28) выберем следующие:

на поверхности Σ₂ известны λ₁, λ₂ и λ₃,

на поверхности сопла σ λ₂ = 0.

В качестве условия оптимальности поверхности сопла σ будем рассматривать условие (26). Таким образом, мы выписали корректную краевую задачу для уравнений (22)–(26).

3. АЛГОРИТМ РАСЧЕТА ОПТИМАЛЬНОЙ ФОРМЫ СВЕРХЗВУКОВОЙ ЧАСТИ ПРОСТРАНСТВЕННОГО СОПЛА

В [9] был разработан численный метод пространственных характеристик для расчета сверхзвуковых течений газа. Входными параметрами задачи являлись числа M – количество меридиональных плоскостей φ = const на отрезке [0, π], N – количество плоскостей x = const, секущих заданную начальную осесимметричную характеристическую поверхность Σ₀, и P – количество характеристических поверхностей в пространственном течении разрежения Прандтля–Майера. Линии пересечения характеристических поверхностей противоположных семейств будем называть венками. Расчетная сетка строилась пересечением венков и меридиональных плоскостей в процессе решения задачи.

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

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

Опишем теперь общую последовательность действий, в соответствии с которой мы будем решать задачу профилирования. Длина сопла L предполагается фиксированной. В плоскости x = L задан гладкий контур γ₂. Оптимальная стенка сопла σ искалась в два этапа, на каждом из которых был организован итерационный процесс. На первом этапе проводился прямой расчет сопла заданной формы r = r(x, φ), определялось распределение узловых точек и находились скорость V, функция тока ψ в узловых точках области τ и функция тока χ в фиксированных точках стенки сопла. Для того чтобы рассматриваемая сетка целиком заполняла область τ, т.е. для того, чтобы замыкающая характеристическая поверхность Σ₂, найденная в результате прямого расчета, проходила через заданный выходной контур γ₂, был организован итерационный процесс, в результате которого на поверхности Σ₁ находился венок γ₀, с которого замыкающая характеристическая поверхность Σ₂ приходила на контур γ₂. Остальные венки на Σ₁ сдвигались пропорционально. При этом количество венков на Σ₁ сохранялось. Газодинамические функции в сдвинутых венках находились с помощью квадратичной интерполяции по данным на исходной поверхности Σ₁.

На втором этапе решалась сопряженная задача. Сначала решалась задача Коши для уравнения (23) с начальными данными (25), и находилась функция λ₁ на замыкающей характеристической поверхности Σ₂. Затем на этой поверхности из (29) определялась функция λ₂. Далее решалась задача Дарбу для системы (28) с известными значениями λ₁, λ₂ и λ₃ на Σ₂ и условием (30) на σ, и определялись значения искомых функций λ₁, λ₂ и λ₃ в узлах найденной ранее сетки в области τ. Как в задаче Коши, так и в задаче Дарбу решение искалось с помощью трех пересчетов. Вычисленные в узловых точках стенки значения λ₁ подставлялись в функцию невязки $E(x,~\varphi ,~{{\lambda }_{1}})$ условия оптимальности (26). Теперь ординаты новых узловых точек стенки в каждой меридиональной плоскости ${{\varphi }_{*}}$ вычислялись по формуле

Здесь $F(x,~\varphi )$ – некоторая заданная функция, удачный выбор которой позволял уменьшить общее число итераций, необходимых для нахождения оптимальной стенки.

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

4. ПРИМЕРЫ РАСЧЕТОВ

Во всех приведенных расчетах $\unicode{230} = 1.4$, ${{p}_{0}} = 0$.

Были проведены два расчета сверхзвуковой части сопла длиной L = 7.6 с круглым критическим сечением r = 1 с выходным контуром в виде прямоугольного сектора круга радиуса R = 6.8974 со скругленными по радиусам ${{r}_{1}}$ и ${{r}_{2}}$ углами. На фиг. 2 показан общий вид выходного контура ${{r}_{L}}(\varphi )$. При использовании входных параметров M × N × P = 31 × 61 × 19$~$относительная погрешность вычисления тяги не превышала 0.1%.

В первом расчете радиус скругления углов был одинаков и равнялся ${{r}_{1}} = {{r}_{2}} = 2.5570$. В качестве начального приближения бралась форма, рассчитанная прямым методом оптимизации (см. [1]). Осесимметричная поверхность Σ₀ бралась из таблиц [11]. Функция $F(x,~\varphi )$ из (31) задавалась в виде $F(x,~\varphi ) = 1 + \beta x$, где параметр β в процессе итераций менялся от 10 до 0.

Тяга оптимальной сверхзвуковой части сопла в направлении единичного вектора ${\mathbf{e}} = - {{{\mathbf{e}}}_{x}}$ составила ${{T}_{x}} = 1.549$, а в поперечном направлении ${{T}_{y}} = - 0.079$. Тяга же оптимальной осесимметричной сверхзвуковой части сопла с выходным контуром, вписанным в данный (${{r}_{{{\text{вп}}}}} = 2.8570$), составила ${{T}_{x}} = 1.522$. Таким образом, абсолютная величина тяги оптимальной сверхзвуковой части пространственного сопла (T = 1.551) превышает тягу закритической части вписанного осесимметричного сопла на 1.9%.

Во втором расчете радиусы скруглений брались разные: ${{r}_{1}} = 2.3$, ${{r}_{2}} = 2.0$. Тяга в продольном направлении составила ${{T}_{x}} = 1.555$, а в поперечном направлении $~{{T}_{y}} = - 0.073$. Таким образом, абсолютная величина тяги этой оптимальной сверхзвуковой части сопла (T = 1.557) превысилa тягу вписанного сопла на 2.2%.

Был проведен расчет сверхзвуковой части сопла длиной L = 14.96 с круглым критическим сечением r = 1 с выходным контуром в виде 120° – сектора круга радиуса R = 12.0711 со скругленными по радиусам ${{r}_{1}} = 5.4$ и ${{r}_{2}} = 4.5$ углами. На фиг. 3 показан общий вид выходного контура ${{r}_{L}}(\varphi )$. Тяга оптимальной сверхзвуковой части сопла составила Т = 1.864. Для сравнения (с использованием той же сетки) был проведен расчет оптимального осесимметричного сопла с ${{r}_{L}}(\varphi ) \equiv 5.6022$, вписанного в рассчитанное пространственное сопло. Его тяга составила $T = 1.855$. Таким образом, величина тяги оптимальной пространственной сверхзвуковой части сопла превысила тягу вписанного осесимметричного сопла на 0.5%.

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

Отметим, что в случае вывода спутника на околоземную орбиту увеличение тяги ракетного сопла на 0.3% позволяет вывести на орбиту на 1.3% больше полезного груза, или поднять высоту орбиты на 10% (см. [12]).