CC BY
28
УЧЕНЫЕ ЗАПИСКИ Ц АГ И Том XXI 19 90
№ 6
УДК532.542 + 532.525.011.5
ОБ ИТЕРАЦИОННЫХ АЛГОРИТМАХ РЕШЕНИЯ ЗАДАЧ ПРОФИЛИРОВАНИЯ ПРОСТРАНСТВЕННЫХ
КАНАЛОВ
Ю. В. Куриленко, М. П. Левин
Рассматриваются задачи профилирования сверхзвуковых пространственных каналов, реализующих на выходе равномерный поток или заданное распределение газодинамических функций на центральном теле. Эти некорректные задачи сводятся к решению ряда последовательных задач обтекания параметризованной стенки. Параметры, определяющие форму стенки, определяются из решения систем нелинейных уравнений. Рассмотрены два итерационных алгоритма решения таких систем. Первый алгоритм основан на использовании метода релаксации, второй — на использовании метода Бройдена. Приведены результаты сопоставления этих алгоритмов на примере решения двух задач профилирования пространственных каналов.
1. Рассмотрим две задачи профилирования сверхзвуковых пространственных каналов, исследовавшихся в работах [1, 2]. Первая задача посвящена профилированию сверхзвуковых частей пространственных сопл, реализующих на выходе равномерный поток с заданным числом Маха за замыкающей характеристической поверхностью пространственного вида.
В качестве замыкающей характеристической поверхности F+ рассматривались поверхности следующего вида. Пусть ось х совпадает с направлением равномерного потока. В плоскости х=х0—const рассмотрим правильный звездообразный л-угольник. Из внутренних углов этого многоугольника выпустим п полубесконечных конусов с полууглом раствора aK=arcsin (1/М). Оси этих конусов параллельны оси х и направлены в сторону, противоположную оси х. Боковая поверхность пирамидообразного тела, получающегося в результате пересечения конусов, представляет собой первую часть замыкающей характеристической поверхности. Вторая часть является частью боковой поверхности пирамиды с плоскими гранями, проходящими через стороны звездообразного я-угольника и касающимися конуса с полууглом раствора ак, вершина которого совпадает с вершиной пирамиды. Подробно обстоятельства, диктующие необходимость такого выбора замыкающей характеристической поверхности, обсуждаются в работе [3]. Пример формы замыкающей характеристической поверхности изображен на рис. 1,а. Начальные данные задавались на осесимметричной характеристической поверхности F- из таблиц [4]. Эта поверхность стыковалась с F+ в точке К, лежащей на оси х. Пусть в начальном сечении х=0 контур сопла Г- является окружностью единичного радиуса с центром в точке x=i/=z=0.
Задача профилирования 1 формулируется следующим образом. В сверхзвуковом потоке (| V | /а> 1), удовлетворяющем уравнениям
div (р V) — 0, rot V = О
и начальным данным на характеристических поверхностях F- и F+, требуется определить поверхность тока F, проходящую через заданный контур Г- (рис. 1, б). Здесь V — вектор скорости, а плотность газа р и скорость звука а являются известными функциями модуля вектора скорости.
Вторая задача посвящена профилированию сверхзвуковой части пространственного канала, реализующего на заданном центральном теле F0 требуемое распределение газодинамических функций (рис. 1,е). Начальные данные задаются на части характеристической поверхности F~, ограниченной контурами Г- и Г. Здесь Г — линия пересечения поверхности F~ с поверхностью центрального тела F0. Задача профилирования 2 формулируется следующим образом. В сверхзвуковом потоке, удовлетворяющем уравнениям
div (р И) = о > V X го* V + Т grad S = 0 , начальным данным на характеристической поверхности F- и условию
М=/(лг,Ф) (1)
на поверхности центрального тела F0, требуется определить поверхность тока F, проходящую через заданный контур Г-. Здесь плотность газа р и абсолютная температура Т являются известными функциями энтропии S и модуля вектора скорости, f — заданная функция, М = [К|/а — число Маха.
Поверхности F+, F~ и F0 в задачах 1 и 2 предполагаются кусочно-гладкими, а функция f—^непрерывной. Гладкость выше перечисленных поверхностей может нарушаться только в меридиональных плоскостях, являющихся плоскостями симметрии.
2. В работах [5—7] показано, что численное решение сформулированных задач профилирования может быть получено в процессе решения серии элементарных задач обтекания параметризованной поверхности тока F.
Рассмотрим вопрос о параметрическом задании поверхности тока F в цилиндрической системе координат х, г, ср и численном решении серии задач обтекания.
Обозначим через Q число плоскостей симметрии в рассматриваемом потоке. Для получения аппроксимирующей системы соотношений /-го приближения введем /+1 меридиональную плоскость
¥¿=2* (k -l)(JQ)-\ 6=1, 2........... J+ 1
на отрезке 0«p<2jr/Q.
Форму поверхности тока Р в каждой меридиональной плоскости на отрезке
Х{] зададим в виде полинома Эрмита третьей степени
Ч = (&11Щ-1 5 О - £)3 - (*|/Л)| (1 - 6) + «а (3 - 26) ,
Ч = (г - /■,_,) (о - 1, 6 «= (Г - Х,_!) (ДГ| - *,_!)- 1.
Величины *<_1, /Ч_1 (|3'П/<5|) г_1, ¿¡, (дт]/д|){ считаются известными. Величины Г*
определяются из решения задачи обтекания участка поверхности тока £ при заданных значениях управляющих параметров (<?т]/д|); методом пространственных характеристик. При решении задачи 1 использовался вариант метода пространственных характеристик, изложенный в работе [8], а при решении задачи 2 — вариант, изложенный в работе [9]. Смешанная задача обтекания решается вплоть до замыкающей характеристической поверхности 5+'в задаче 1 и до поверхности центрального тела £° в задаче 2.
Выбор управляющих параметров (<5т]/д£){, определяющих форму поверхности тока Р, в задаче 1 осуществлялся из удовлетворения конечно-разностной аппроксимации условия совместности вдоль характеристической поверхности второго семейства (приходящей с поверхности тока Р). Выполнение условия совместности было необходимо в точках пересечения этой характеристической поверхности с замыкающей характеристической поверхностью £+ (первого семейства). Поскольку поверхность /•"+ отделяет область равномерного потока, то рассматриваемое условие совместности в дифференциальной форме представлено в виде [6]:
В (— 2> п—. Я—) I (¿и!<1л:) ех + (¿и/«Гдт) е г+ (йт\йх) е?| = 0 ,
2 = (¿*/¿9) ех + (¿г/<*р) ег +г (2)
В (а, Ь, с) = а X Ъ — М (М а Ь) + Ь (М с а).
Здесь я_— единичный вектор, направленный по нормали к характеристической поверхности второго семейства, М= У/а — вектор Маха, полные производные по х вычисляются вдоль линий пересечения характеристической поверхности второго семейства с меридиональными плоскостями, полные производные ПО ф вычисляются вдоль линии пересечения характеристической поверхности второго семейства с поверхностью £+, ех, ег и е,_ — орты цилиндрической системы координат.
В задаче 2 выбор управляющих параметров (<5т]/д|) осуществлялся из удовлетворения условию (1) в точках поверхности центрального тела £°.
3. Обозначим через Егь, (&=1, 2,...,/+1) невязки конечно-разностной аппроксимации условия совместности (2) в различных меридиональных плоскостях для задачи 1 или невязки условия (1) для задачи 2, соответствующие заданным параметрам (<?Т1/<Э\)1к. На практике точные значения невязок Ец1 вычислить невозможно. Погрешность вычисления невязок £¡4 определяется погрешностью численного метода пространственных характеристик, используемого для решения задач обтекания, и, кроме того, для задачи 1 — погрешностью конечно-разностной аппроксимации условия совместности (2). Таким образом, параметры, определяющие форму стенки (<Эт1/с)£)гь, должны удовлетворять нелинейной системе уравнений. Формально эту систему можно записать в следующем виде
Ог,- = 0, (3)
где О — нелинейный оператор, а г/ = { (<?ч/д£г)/й. к — 1, 2,..., 7+1}. Отметим, что конкретный вид оператора О нам неизвестен. Мы только можем вычислять вектор невязки £¡ = 0^, где £г={£;1, Е12.............£/)(/+1)}.
Применение классического метода Ньютона для решения системы (3) невозможно, поскольку нельзя аналитически вычислить якобиан этой системы. Использование конечно-разностных аппроксимаций для вычисления якобиана, как показывает опыт расчетов, приводит к потере точности вычислений из-за ошибок аппроксимации, входящих в компоненты вектора невязок.
В работе [6] для решения задачи (3) использовался следующий алгоритм (алгоритм А)-.
г\+х = г[ + В} Е{.
Здесь У — номер итерации. В* — диагональная матрица с элементами ь£, вычисляемыми по формулам
bi I ifi > hi|<c’
* C sign (ц^) , I Tfft I > C *
Л = - {(Л|/«)І - (дпЩІҐ] / И* - ЕІҐ) •
Постоянная С подбирается экспериментально так, чтобы обеспечить сходимость итерационного процесса.
Для решения системы (3), наряду с алгоритмом А, используем метод Бройдена [10, 11]. В этом случае итерационный алгоритм (алгоритм 5) решения системы (3) имеет следующий вид: для известных на /-й итерации значений якобиана системы (3) А> и вектора невязки Е[ решаем систему линейных уравнений
На первой итерации якобиан задается в виде диагональной матрицы с одинаковыми ненулевыми элементами. В процессе итераций матрица якобиана становится полностью заполненной и несимметричной, что следует иметь ввиду при решении системы линейных уравнений на шаге 1 в алгоритме Б.
4. Рассмотрим результаты применения алгоритмов Л и £ при решении задач профилирования 1 й 2.
При численном решении задачи 1 рассматривалось течение с Q—8 плоскостями симметрии, реализующее равномерный поток с числом Маха М=4,273. Расчеты проводились на характеристической сетке при 1=6. Число опорных точек N, определяющих искомую поверхность тока F в каждой меридиональной плоскости, было равным 58. При этом решалось 56 элементарных задач обтекания. Форма полученной в процессе расчета поверхности тока F изображена на рис. 2 и 3. На рис 2 представлены образующие поверхности тока F в трех меридиональных плоскостях, а на рис. 3 поперечные сечения этой поверхности плоскостями х=const. Поскольку замыкающий контур поверхности тока F не лежит в плоскости x=const, при изображении поперечных сечений поверхность F дополнялась воображаемой цилиндрической поверхностью с образующей параллельной оси х. На рис. 4 приведен график зависимости числа итераций / в алгоритмах А я Б при решении элементарных задач обтекания в зависимости от порядкового номера задачи i (/=2, 3,..., N—1), совпадающего с номером опорной точки на искомой поверхности тока F.
Из представленных результатов видно, что число итераций при решении элементарных задач обтекания с помощью алгоритма Б меньше, чем при использовании алгоритма А. Однако число арифметических операций на каждой итерации в алгоритме Б больше, чем в алгоритме А, поскольку приходится решать систему (/+1) линейных уравнений с полностью заполненной матрицей коэффициентов. Таким образом, например, при использовании QR-алгоритма из работы [12] в алгоритме Б требуется порядка 2(/+1)3/3 дополнительных операций умножения и деления. Однако при рассматриваемых значениях J число дополнительных операций умножения и деления намного меньше числа операций, требующихся для вычисления вектора невязки Е. В результате полное время решения задачи профилирования 1 на рассматриваемой сетке с помощью алгоритма Б на ЭВМ БЭОМ-6 оказалось заметно меньше (1 час 15 минут) времени решения с помощью алгоритма А (1 час 40 минут).
В задаче 2 в качестве центрального тела F0 рассматривался круговой конус с полууглом раствора 15°. Распределение числа Маха на конусе и начальные данные на характеристической поверхности F~ задавались из решения задачи обтекания конуса под углом атаки 10° сверхзвуковым потоком газа с показателем адиабаты х=1,4 и значением числа Маха в невозмущенном потоке М,»=3. Контур Г- задавался аналитически
Здесь Я (<р) — координата г точки, лежащей на ударной волне в меридиональной плоскости ф=сопз1 при х=1. В расчетах значения /?(ф) брались из работы [13].
Расчеты проводились на характеристической сетке при /=16. Число опорных точек, определяющих искомую поверхность тока 7\ было равным 21.
*(?)=!, г (?) = 0,5 [tg 15° + /? (?)] .
«Ученые записки» № 6
Решение рассматриваемой задачи с помощью алгоритма А получить не удалось яз-за отсутствия сходимости итерационного процесса. Применение алгоритма Б оказалось успешным. Для достижения необходимой точности требовалось от 5 до 10 итераций при решении элементарных задач обтекания. На рис. 5 представлены образующие найденной поверхности тока F в трех меридиональных плоскостях (<р=0, <р = хс/2 и ср=я). Образующая поверхности тока F, лежащая в плоскости •Ф=0, соответствует наветренной стороне конуса, а образующая, лежащая в плоскости <ф=я соответствует подветренной стороне конуса.
ЛИТЕРАТУРА
1. Борисов В. М., Левин М. П. Расчет и оптимизация сверхзвуковых частей пространственных сопл с косым срезом. — М.: ВЦ АН СССР, 1987.
2. Куриленко Ю. В., Михайлов И. Е. Профилирование сверхзвуковых пространственных каналов с центральным телом. — М.:
ВЦ АН СССР, 1989.
3. Б о р и с о в В. М. О структуре решений при оптимизации формы сверхзвуковых пространственных каналов.—М.: ВЦ АН СССР, 1987.
4. Кацкова О. Н., Шмыглевский Ю. Д. Таблицы параметров осесимметричного течения свободно расширяющегося газа с плоской переходной поверхностью. — М.: Изд-во АН СССР, 1962.
5. JI е в и н М. П. Об оптимальных поверхностях тока в сверхзвуковых пространственных течениях. — Ж. вычисл. матем. и матем. физ. —
1982. — Т. 22, № 4.
6. J1 е в и н М. П. Метод контрольной поверхности в пространственных вариационных задачах газовой динамики. — М.: ВЦ АН СССР, 1983.
7. Михайлов И. Е. О численном решении пространственного аналога задачи Гурса. — Ж- вычисл. матем. и матем. физ.—'1983.—
Т. 23, № 5.
8. Л е в и н М. П. Схема метода пространственных характеристик с использованием бихарактеристических соотношений. — М.: ВЦ АН СССР,
>1985.
9. Борисов В. М., Куриленко Ю. В., Михайлов И. Е., Николаевская Е. Л. Метод характеристик для расчета вихревых сверхзвуковых пространственных установившихся течений. — М.:
ВЦ АН СССР, 1988.
10. В г о у d е п С. G. A class of methods for solving nonlinear si-multaneous équations. Math, of Comput. — 1965. N 19.
11. Grzegorsky S. M. Orthogonal projections on convex sets for Newton—like methods.—SIAM J. Numerical Anal. — 1985. V. 22. N 6.
12. Дэннис Д ж., Шнабель P. Численные методы безусловной «оптимизации и решения нелинейных уравнений. — М.: Мир, 1989.
13. Б а б е н к о К. И., В о ск р е с е н ски й Г. П., Любимов А. Н., Русанов В. В. Пространственное обтекание гладких тел идеальным газом. — М.: Наука, 1964.
Рукопись поступила 13/Х 1989 г.
