Определение формы профиля по заданной хордовой диаграмме чисел Маха в трансзвуковом потоке Текст научной статьи по специальности «Математика»

Т о м XV

УЧЕНЫЕ ЗАПИСКИ Ц А Г И

19 84

№ 4

УДК 533.6.011.3:629.7.025.73

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

А. А. Шагаев

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

Один из возможных подходов к проектированию крыловых профилей с желаемыми аэродинамическими характеристиками основан на решении обратной задачи с заданной хордовой диаграммой, т. е. на определении формы профиля по заданному вдоль хорды распределению параметров потока на его поверхности. В работах [1—3] описаны итерационные методы решения этой задачи для несжимаемого течения, а в работах [4—6] — разностные методы, применимые в случае течений с большими дозвуковыми и трансзвуковыми скоростями. Следует отметить, что реализация метода [4] возможна только в режиме диалога с ЭВМ. В методе [5] решение ищется в классе разомкнутых самопересекающихся контуров и физический смысл получаемых при этом решений не всегда ясен, на что указано, в частности, в работе [6], где обратная задача рассмотрена в приближении теории малых возмущений.

Метод последовательных приближений |[1] свободен от указанных недостатков и позволяет находить приемлемые решения обратной задачи даже в случае некорректно заданного распределения параметров потока- Обобщение этого метода на случай течений газа С. А. Чаплыгина |[7] дано в работе ¡[8], а в настоящей работе излагается его дальнейшее развитие на случай дозвуковых и трансзвуковых течений совершенного газа.

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

4 Ф,ФуФ,у = 0, (1)

х+1 * *

где к — отношение удельных теплоемкостей, и удовлетворяет на поверхности профиля у = !±(х), хе [0, 1] условию непротекания

=0. (2)

дп у=/±(л:)

Здесь функция /+ {х) соответствует верхней, /"(*)—нижней поверхности профиля, а д/дп — производная по нормали к его поверхности.

Распределение коэффициента скорости % на поверхности профиля выражается через распределение потенциала следующим образом:

— -х, (3)

д1 у=/±(*)

где д/дI — производная по касательной к поверхности профиля.

При х2+ г/2—>-оо производные потенциала ограничены и связаны с коэффициентом скорости набегающего потока Хос и углом атаки а (углом наклона вектора скорости набегающего потока к оси х) соотношениями

= + Фу, (4а)

Ф

(46)

"X

Обратная задача с заданной хордовой диаграммой состоит в определении функции Ф(х, у) и области, в которой эта функция удовлетворяет уравнению для потенциала (1), т. е. требуется найти потенциал Ф и функции ^±(х), описывающие форму профиля. Для этого на поверхности искомого профиля задается на одно граничное условие больше, чем в случае прямой задачи, а именно условие непротекания

(2) и распределение коэффициента скорости вдоль хорды {условие

(3)1.

В плоскости (х, -ф), где — функция тока, задача (1) — (3) сводится к задаче Дирихле для некоторого квазилинейного уравнения (см. монографию [9]) относительно функции %, которая определяется соотношением — где р — давление, а V — составляющая

вектора скорости вдоль оси у.

В работе |[10] была показана однозначная разрешимость задачи (1) — (3) в классе разомкнутых и быть может самопересекающихся контуров для равномерно дозвуковых течений без точек торможения. При этом параметры набегающего потока и а не задаются, а определяются из решения обратной задачи. По-видимому, аналогичное утверждение справедливо и в более интересном для практики случае течений около крылового профиля с затупленной передней кромкой, на что указывают результаты исследования обратной задачи с заданной хордовой диаграммой в рамках дозвуковой линейной теории, а также обратной задачи, когда носителем данных является длина дуги поверхности профиля [2].

При проектировании крыловых профилей требуется, чтобы желаемая хордовая диаграмма реализовалась при заданном значении Яоо- С математической точки зрения это означает, что задачу (1) — (3) требуется решать с дополнительным условием (4а). Для существова-

ния такого решения заданное распределение К(х) должно удовлетворять некоторому условию, которое вместе с условиями замкнутости контура относится к условиям разрешимости [2]. Однако их выполнение не гарантирует еще однолистности получаемого контура.

Таким образом, принципиальное отличие обратной задачи (1)— (4а) от прямой заключается в некорректности обратной задачи в том смысле, что для произвольно заданного распределения чисел М, М3(х) или коэффициента скорости %3(х) вдоль поверхности искомого профиля ее решение может либо не существовать, либо соответствовать физически нереальному контуру. При этом проверка условий разрешимости эквивалентна решению задачи (1) — (3). Поэтому для решения обратной задачи с заданной хордовой диаграммой целесообразно применять такие итерационные методы, в которых каждое итерационное приближение удовлетворяет условию (4а) и описывает течение около замкнутого однолистного профиля.

В случае задачи с заданной хордовой диаграммой скорости У3(х), где скорость V отнесена к ее значению в набегающем потоке, условие (4а) записывается в виде V х2+у2^оо, применимом при Моо = 0.

Естественно потребовать, чтобы рассматриваемый итерационный процесс сходился, если заданная хордовая диаграмма реализуется на крыловом профиле при заданных параметрах набегающего потока. В случае несжимаемого течения таким свойством обладает метод работы [1]. Кроме того, расчеты показывают, что в этом методе наиболее быстро по итерациям заданное распределение скорости устанавливается на основной части хорды. Поэтому с его помощью можно получать приемлемые для проектировщика решения даже в случае произвольно заданного распределения скорости за счет того, что распределение скорости на профиле, соответствующее текущей итерации, отличается от заданного лишь в окрестности передней и задней кромок. В настоящей работе аналогичный итерационный процесс решения обратной задачи предлагается для дозвуковых и трансзвуковых течений.

2. Итерационный метод. Итерационный процесс решения обратной задачи с заданной хордовой диаграммой заключается в определении некоторой постоянной ю0 и функции ан+гЧь аналитической внутри круга единичного радиуса вспомогательной плоскости £ и удовлетворяющей условиям

М0) = 0, Т1(1) = 2я. (5)

Для определения формы профиля по функции <В1(£)-И'т;1(£) и постоянной со0 следует при помощи отображающей функции

= е-.<ч+*ч<о (6)

конформно отобразить единичный круг |С|-<1 на внешность близкой к окружности кривой в плоскости г1—ею+к, где (со, %) — криволинейные координаты — оо < № < сю, О т -< 2к. Координаты профиля в физической плоскости г = х -¡- 1у определяются по криволинейным координатам кривой ш(т), на которую отображается окружность 0 <>-<2«, при помощи простых соотношений

1+С08[т(е)] .. 8Ь[ш(е)]81пК0]

2 ' 2 ' (7)

где т(а) = т!(е) — е, а ш (в) = со1 (е) + со0 в соответствии с выражением (6).

2—«Ученые записки» № 4

17

Выражения (6) и (7) устанавливают взаимно однозначное соответствие между окружностью |£| = 1 и произвольным замкнутым контуром у=/= (х) с единичной хордой, расположенной на отрезке оси х при 0сл;с1, если функции ¡±(х) удовлетворяют условию

lim ^jfL-Hm-^

х->1 VI— х х-1 VI— X

Отметим, что это условие выполняется для рассматриваемых в аэродинамике крыловых профилей; имеющих угловую или затупленную (с конечным радиусом кривизны) заднюю кромку. Если известна форма такого профиля, то соответствие между точками его поверхности к окружности |£| = 1 определяется в результате конформного отображения единичного круга на внешность кривой (о = со(т)

ш (т) = In

У*™ , -i/i + i

Sin X V Sin2 X

(8)

где T=2arccos(}/x) на верхней поверхности и х = 2 [тс — arccos (Ух)) на нижней поверхности профиля.

Это конформное отображение, а следовательно, и постоянная <о0 вместе с удовлетворяющей условиям (5) аналитической функцией u>j (С) + iix (С) определяются единственным образом в соответствии с теоремой Римана [11J, если потребовать, чтобы задняя кромка профиля т = 0,2я отображалась в точку С = 1.

Для определения формы профиля по распределению скорости V(х) = sgn К — т) У (т)> 1/(хо) = 0. х0^(0,2те) на его поверхности достаточно знать функцию

<■>„ • в0~£ • е е 0 sin —-— sin —

---í-JL (9)

v V(z)$(V, Moo) dx v '

где £q — угловая координата передней точки торможения в плоскости С, а функция {3(1/, Mo,) выбирается ниже из условия наиболее быстрой сходимости итераций при Моо > 0, ,6(1/, 0)=1.

Постоянные е0, <в0 и зависимость т(е), удовлетворяющую условиям T(0) = 2it, т (г0) = т0, т (2iu) = 0, можно определить путем решения дифференциального уравнения

F (х) ~~~~ — еЮ°sin sin — , (10)

W de 2 2 V '

которое получается из равенства (9) при F(z) = D(x) l/(t){3(V, v). Затем по формуле (6) можно найти значения функции -^(С) на окружности í = eh и, решая задачу Шварца [12] для круга |С|<1, определить аналитическую функцию cu^O + í-c,^), а следовательно, и форму профиля.

Таким же образом, как и в работе [1], можно показать, что с помощью описанной выше процедуры для любой непрерывной функции F(i), которая один раз меняет знак на интервале (0,2т:) и удовлетворяет условию F(0) — F(2k), однозначно определяется замкнутый профиль. На этом и основан итерационный процесс решения обратной задачи. По некоторой функции D(x) и заданному распределению скорости 1/3(*) можно вычислить функцию F(i) —

= D (x) 1/3 (x) fJ ( 1/3, M«,) и найти соответствующий ей контур. Если функции D(x) и К3(х) отвечают одному и тому же профилю, то таким образом восстановится его форма. В противном случае получится новый контур и определится его отображение (6), (7) на окружность С = е'е. Для замыкания процесса итераций нужно рассчитать обтекание этого контура и найти новую функцию £>(х) по формуле (9).

При Моо = 0 описанный выше итерационный процесс переходит в процесс работы [1], если соотношения (7) заменить преобразованием Жуковского

Тогда отображение (6), (11) оказывается конформным и функцию D (т) можно вычислять по формуле

D (х) = -у /(sh2 <0 + sin2 x) (1 +0,2), (12)

так как в этом случае выражение (9) совпадает с известной формулой Теодорсена—Серебрийского |[13] для распределения скорости на поверхности профиля в несжимаемом течении. Сходимость итерационного процесса 1[1] обусловлена тем, что форма получающихся на итерациях контуров определяется главным образом заданным распределением скорости, поскольку зависимости Di(-r) для разных профилей различаются между собой членами О (со2), где со — величина порядка толщины профиля. Определенные соотношениями (6), (11) и (6), (7) отображающие функции также отличаются малыми величинами 0(ю2), поэтому следует ожидать, что при Моо = 0 итерационный процесс настоящей работы сходится тогда, когда сходится итерационный процесс работы |[1]. В то же время применение отображающей функции (6), (7) позволяет зафиксировать по итерациям зависимость V3(t) и положение передней кромки профиля в плоскости Zi, что существенно упрощает логику вычислительного алгоритма.

Описанный процесс решения обратной задачи разделяется на две практически не зависящие друг от друга части. В первой из них определяется форма контура, а во второй рассчитывается распределение скорости по этому контуру, причем от способа решения задачи обтекания итерационный процесс не зависит. В случае малых чисел Моо можно, как и в несжимаемом течении, положить Р (V, Moo) = 1 и рассчитывать обтекание профиля сжимаемой жидкостью. Тогда в силу непрерывной зависимости решения от параметра Моо получаемые распределения V (т) будут мало отличаться от соответствующих зависимостей в несжимаемом потоке, что гарантирует сходимость итерационного процесса и в некотором диапазоне чисел Моо, отличных от нуля.

Рис. 1

Вычислительный алгоритм решения обратной задачи изображен на рис. 1. Для начала решения задается нулевое приближение £>о(т). Затем по функции Fn (т) = £>n-i (т) V3(t) Р^.М«,) путем решения дифференциального уравнения (10) определяются зависимость т(е) и постоянная и0, как показано в работе [8], а по формуле обращения Гильберта [12] находятся краевые значения функции k>i(£) на окружности 1^ = еи. Координаты профиля вычисляются по формулам (7). Перед решением задачи обтекания проверяется условие однолистности контура со(т)+о)(2я—т)>0. Если это неравенство нарушается, то постоянная cûo увеличивается на минимальную константу Дсоо, необходимую для его выполнения. Полученный таким образом однолистный контур соответствует прежним зависимостям coi(e) и т(е).

Обтекание профиля рассчитывается с заданным значением коэффициента подъемной силы су, величина которого определяется по известному из постановки задачи значению коэффициента нормальной силы Cjv и углу атаки профиля с предыдущей итерации an-i(a0 = O). Вычислением функции Dn(т) по формуле (9) итерационный процесс замыкается. Если в качестве начального приближения задается какой-либо исходный контур {х), то для начала решения обратной за-

дачи требуется найти его отображение (6), (7) на окружность |£| = 1 путем конформного отображения внешности кривой со(т) (8) в плоскости Zi на единичный круг в плоскости Ç, для чего применяется алгоритм, основанный на методе работы [13].

Описанный алгоритм был реализован в виде программы для ЭВМ БЭСМ-6 на автокоде БЕМШ. Разработанная программа имеет две версии, которые различаются алгоритмами решения задачи обтекания профиля. В первой версии используется алгоритм, основанный на быстром методе расчета дозвуковых сжимаемых течений около профиля в приближении газа Чаплыгина [14], а во второй версии — алгоритм решения полного уравнения потенциала (1) релаксационным методом конечных разностей [15].

Проведенные на ЭВМ расчеты показали, что при р=1 итерационный процесс решения обратной задачи сходится в случае несжимаемой жидкости за несколько итераций, а с ростом числа Мао его сходимость замедляется. Скорость сходимости сохранится, если в итерационном процессе использовать такую функцию Р(К, Мое), при которой процесс определения формы профиля по распределению V3{x) в сжимаемой жидкости будет совпадать с процессом определения формы этого же профиля по распределению VH.3(х) =р(У3, Мю) V3(x) в несжимаемой жидкости. Для этого на поверхности профиля должно выполняться равенство VH=p(V, М^) V. В общем случае это равенство не имеет места, однако в приближении тонкого профиля можно установить его аналог

где р — плотность, отнесенная к ее значению в набегающем потоке.

Таким образом, для дозвуковых сжимаемых течений оптимальной является следующая функция р:

В случае течений с местной сверхзвуковой зоной для функции р сохранялся лишь качественный характер зависимости (13)

(13)

р=рт.

(14)

Для того, чтобы при таком выборе функции р итерационный процесс решения обратной задачи в сжимаемой жидкости был эквивалентен в рамках дозвуковой теории тонкого профиля итерационному процессу определения аффинно-подобного контура по распределению скорости р? (V3, Moo) V3(t) в несжимаемой жидкости, отображение (7) следует изменить:

1 + COST У* ~~~ Мте sh [ы (е)] sin [т: (е)]

Расчетные исследования показали, что в случае трансзвуковых течений со скачками уплотнения для устойчивости итерационного процесса решения обратной задачи нужно ввести релаксацию функции D:

■V-

~2— sin ~2~

VnMP(V«. MJ

djn

dt

а постоянную у в выражениях (14) и (15) определять при помощи соотношения

__1

Т_(Мтах + Д)з '

где Мтах — максимальное значение числа М в заданном распределении М3(х), а Л «0,05.

Для течений с местными числами М меньше 1,2 такой итерационный процесс обычно сходится за 3—6 итераций в зависимости от начального приближения. При этом затраты машинного времени на одну итерацию составляет менее 10 с без учета времени, необходимого для решения задачи обтекания профиля.

3. Результаты расчетов. Точность вычислительного алгоритма проверялась путем решения обратной задачи для корректно заданных распределений Мз(х), которые были получены в результате расчета обтекания известных профилей релаксационным методом [15] и заведомо удовлетворяют условиям разрешимости. На рис. 2 и 3 приводятся полученные таким образом результаты соответственно для профиля из [16] и сверхкритического профиля, эпюрный чертеж которого имеется в работе [17]. Распределения М(л:) для этих профилей изображены точками на верхних, а сами профили — на нижних графиках рисунков. Там же для сравнения приведены контуры профилей, получившиеся на 4-й итерации обратной задачи, и рассчитанные для них распределения М(х). При этом в обоих случаях в качестве начального приближения задавалась функция D(т), вычисленная по формуле (12) для руля Жуковского с относительной толщиной 11,98%.

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

На рис. 4 изображены результаты применения настоящего метода для нахождения формы профиля, соответствующего «произвольно заданной» хордовой диаграмме М(х). Хорошо видно, что в результате

распределению М(х) при М<» = 0,75; су = 0,4

Рис. 3. ВосстанФвление сверхкритического профиля по известному распределению М(л:) при М„о=0,8; cv=0,5

Рис. 4

решения обратной задачи получился профиль с близким к заданному «полочным» распределением М.(х) и слабым скачком уплотнения на рассматриваемом режиме течения.

В заключение автор благодарит Ю. Б. Лифшица за руководство настоящей работой.

ЛИТЕРАТУРА

1. Шурыгин В. М. Определение контура профиля по заданному распределению давлений. —Труды ЦАГИ, 1945, вып. 576.

2. Т у м а ш е в Г. Г., H у ж и н М. Т. Обратные краевые задачи и их приложения. — Изд. Казанского университета, 1965.

3. Павловец Г. А., Самознаев Н. Д. Численный метод построения контура крылового профиля по заданному распределению скоростей на его поверхности. — Труды ЦАГИ, 1970. вып. 1271.

4. Tranen T. L. A rapid computer aided transonic airfoil design method. — AIAA Paper N 74—501, 1974.

5. С a r 1 s о n L. A. Transonic airfoil design using. Cartesian coordinates. — NASA CR-2575, 1976.

6. Shankar V., M a 1 muth N. D. and Cole J. D. Computational transonic airfoil design in free air and a wind tunnel. — AIAA Paper N 78—103, 1978.

7. Чаплыгин С. А. О газовых струях. Собр. соч. т. II. — М.: Гостехиздат, 1948.

8. Ш а г а е в А. А. Построение контура профиля по заданному распределению давлений в сжимаемом потоке газа. — Труды ЦАГИ, 1978, вып. 1925.

9. К о ч и н H. Е., К и б е л ь И. А. и Розе Н. В. Теоретическая гидромеханика, ч. И. — М. — Л.;' Гостехиздат, 1948.

10. Монахов В. Н. Об однозначной разрешимости плоских задач газовой динамики со свободными границами. — ДАН СССР, т. 164, № 5, 1965.

11. Евграфов М. А. Аналитические функции. — М.: Наука, 1965.

12. Г ах ов Ф. Д. Краевые задачи. — М.: Физматгиз, 1958.

13. Серебрийский Я. М. Обтекание крыловых профилей произвольной формы. — Труды ЦАГИ, 1944, вып. 553.

14. Ш а г а е в А. А. Приближенный метод расчета дозвуковых сжимаемых течений около профиля. — Ученые записки ЦАГИ, 1976, т. IV, № 5.

15. Лифшиц Ю. Б. К теории трансзвуковых течений около профиля. — Ученые записки ЦАГИ, 1973, т. IV, № 5.

16. Bauer F., Garabedian Р., Korn D. Supercritical wing sections, I, Springer, 1972.

17. Боксер В. Д. Приближенные способы определения начала резкого возрастания сопротивления профиля при околозвуковых скоростях. — Ученые записки ЦАГИ, 1980, т. XI, № 3.

Рукопись поступила 17/XII 1982 г.