Том ХЫН
УЧЕНЫЕ ЗАПИСКИ ЦАГИ 2012
№ 6
УДК 629.735.33.015.3.025.73 533.6.011.35:629.7.025.73
ОБОБЩЕННАЯ ОБРАТНАЯ ЗАДАЧА ДЛЯ ПРОФИЛЯ
А. А. НИКОЛЬСКИЙ
Для общего случая потенциальных течений решена задача о построении контура профиля по заданному на его поверхности распределению давлений и заданной величине продольного момента при нулевой подъемной силе. Для случая обтекания несжимаемой жидкостью получено решение задачи в квадратурах. Для случая обтекания сжимаемой жидкостью, в том числе околозвуковым потоком, построен алгоритм и разработан численный метод решения задачи. Рассмотрен пример решения обобщенной обратной задачи при проектировании вертолетных профилей. Практическая применимость метода подтверждена расчетами на основе СЕБ ЯА^.
Ключевые слова: профиль, обратная задача.
ВВЕДЕНИЕ
Метод решения обобщенной обратной задачи, предложенный в настоящей работе, основан на применении теории функций комплексного переменного. Как известно, обратная задача является, вообще говоря, некорректной. А именно целевое распределение давлений не обязательно соответствует течению с заданной скоростью на бесконечности и/или получаемый в ходе решения контур профиля является замкнутым. Условие сохранения скорости на бесконечности при конформном отображении и два условия замкнутости контура профиля были получены в виде интегральных соотношений в основополагающей работе [1] для несжимаемых потенциальных течений около профиля с нулевой подъемной силой. Далее, как это принято в литературе, эти интегральные соотношения называются условиями разрешимости обратной задачи. В работе [2] получены аналогичные условия для несущих профилей. В работах [3, 4] были предложены методы, позволяющие выполнить эти интегральные ограничения и построить решение обратной задачи. В работах [5, 6] предложены методы решения обратной задачи в рамках линейной теории
дозвуковых течений. Известен также метод работы [7], решающий не- ____________________________
линейную обратную задачу, где условия разрешимости приближенно удовлетворяются в ходе решения с помощью специально подобранных функций. В работе [8] была предложена альтернативная форма интегральных условий разрешимости, а в [9] — метод решения обратной задачи для общего случая нелинейных потенциальных течений, включая течения с внутренними ударными волнами умеренной интенсивности.
Будем называть обобщенной обратной задачей классическую обратную задачу с дополнительным ограничением на величину продольного момента при нулевой подъемной силе.
1. Пусть в физической плоскости г расположен замкнутый контур .
^ ^ „ , _ „ НИКОЛЬСКИМ
заранее не известной формы, обтекаемый потоком идеальной жидко- Александр АлексаНдрОВиЧ
сти. Зададим распределение коэффициента давления ср в зависимости кандидат физико-
р математических наук,
от длины дуги 5, отсчитываемой от задней кромки профиля в направ- начальник сектора ЦАГИ
лении по часовой стрелке. Учитывая, что для потенциальных течений существует взаимно однозначная связь давления и модуля скорости, можно считать заданным и распределение скорости и:
и =
Здесь модуль скорости отнесен к скорости невозмущенного потока, а длина дуги — к полной длине дуги $0 контура профиля. Требуется найти контур профиля, имеющего заданную толщину задней кромки, распределение скорости по которому максимально близко к заданному и при этом имеющий заданную величину продольного момента ет0) при нулевой подъемной силе.
Известно, что решение прямой задачи обтекания профиля потенциальным потоком несжимаемой жидкости находят с помощью конформного отображения внешности контура профиля в плоскости г на внешность окружности в плоскости с'. Для удобства рассмотрения нелинейного случая прибавим к этому преобразование инверсии с = 1/с'. Не ограничивая общности, можно считать, что окружность имеет единичный радиус. Следуя [10], представим производную отображающей функции в виде
где в — угол раствора задней кромки (угол между касательными к верхнему и нижнему контуру) профиля, выраженный в долях п; сп = ап + \Ъп.
При этом действительная и мнимая части логарифма производной отображающей функции могут быть записаны следующим образом:
Здесь ^ — угол наклона касательной к поверхности профиля; 0 — полярный угол в плоскости с.
Для потенциальных течений модуль скорости в плоскости г связан с модулем скорости в расчетной плоскости с следующим соотношением:
Интегрируя это соотношение, получаем зависимость $ (0), находим ^ и определяем координаты контура профиля х (0) и У (0):
(1)
(2)
Как только скорость и (0) известна, решение (2) сводится к квадратурам:
(3)
0
0
0
Осталось выяснить, каким условиям должна удовлетворять функция H = ds для того, что-
d 0
бы решение существовало, получаемый в результате решения профиль имел заданную толщину задней кромки, а его продольный момент при нулевой подъемной силе имел заданную величину.
Интегрируя соотношение (1) и используя теорию вычетов, получим условия сохранения скорости на бесконечности и заданной толщины задней кромки профиля применительно к функции H:
2п 2п 2п
J lnHd0 = 0, J lnHsin(0)d0 = 0, J lnHcos(0)d0 = cn. (5)
0 0 0
Здесь c11 = (1 - s) ——; Ayte — толщина задней кромки, отнесенная к хорде профиля.
2п
Из формулы Чаплыгина — Блазиуса следует дополнительное соотношение для величины продольного момента cm0) при нулевой подъемной силе:
2п
J ln H sin (20) d0 = -Cm0. (6)
0 4
Выписанные интегральные соотношения, как и условия Лайтхилла [1], не позволяют модифицировать целевое распределение скорости до решения задачи, так как для этого требуется
знать зависимость s (0). Однако условия (5), (6) имеют то преимущество, что они применимы
к нелинейному случаю, и далее указан простой способ их удовлетворения.
Для этого, как и в [9], зададим Hm в следующем виде:
Hm = H exp g. (7)
Здесь и далее индекс m означает модифицированную функцию. Будем искать функцию Hm, близкую к заданной функции H в смысле наименьших квадратов, удовлетворяющую интегральным соотношениям (5), (6). Перепишем эти соотношения для функции g:
02 2п
J gd0 = - J ln Hd0 = -a0 2п,
01 0
02 2п
J gcos0d0 = -J lnHcos0d0 = -(a1 -c11 )n,
01 0
02 2n
J g sin 0d0 = -J ln H sin 0d0 = -b, n,
0
2n
ь1п,
0
2n
nHsin20d0 = -f b + -
(8)
J gsin20d0 = -J lnHsin20d0 = -j b2 +Cm0 Jn.
0j 0 4
Нетрудно видеть, что для случая, когда функция H, а с ней и целевая скорость модифицируется по всей поверхности профиля, что соответствует диапазону 0 < 0 < 2п, существует точное решение системы уравнений, имеющее следующий вид:
g = -a0 -(a1 -c11)cos0-b1sin0-jb2 + Jsin20. (9)
0
Это решение отличается от решения для классической обратной задачи [9] наличием последнего слагаемого в правой части соотношения (9). Теперь, подставляя (9) в (7), определяем Нт, а значит и ут. После чего из соотношений (4) находим координаты искомого профиля в параметрическом виде.
В обратной задаче, в отличие от прямой, угол атаки заранее неизвестен, а определяется из
соотношения а = /(у). Здесь у = Г ; Гт;п — циркуляция скорости вдоль части контура про-
Гтш
филя от задней кромки до передней точки торможения вдоль нижней поверхности профиля; Г — полная циркуляция:
еп
2п
Гтш = j и (е) d е, Г= j и (е) d е.
0 0
Для несжимаемой жидкости зависимость скорости от угла атаки известна: и (е) =-2 "sin (е + а) + sin aj, ео = п-2а. Выполнив интегрирование, получим следующее уравнение для определения а:
tg а =
Y
[пу/2 -1 -ау]
Это означает, что для случая несжимаемой жидкости решение обобщенной обратной задачи получено в квадратурах. В нелинейном случае решение и угол атаки получаем на текущей итерации, о чем будет сказано далее.
Наибольший интерес представляет решение с сохранением распределения давлений по верхней поверхности профиля, определяющего уровень максимального коэффициента подъемной силы Су тах. В этом случае модификация выполняется вдоль нижней части контура профиля, и решение не является универсальным. Зададим диапазон модификации 0! <0<02, обо-0 — 0
значим ф =--, тогда 0 < ф< 2п. Как и в работе [9], будем искать функцию g в виде конечного
2п
ряда Фурье:
N
g = ^ dn cos «ф + e« sin «ф.
(10)
«=0
На концах интервала модификации накладываем условия гладкого сопряжения: g (0) = g (2п) = 0, g'(0) = g'(2п) = 0. Решение методом наименьших квадратов сводится к поиску минимума функционала:
02 (02
J g 2dе + А1 g (0) + А,2g'(0) + А,3 J g sin еd0 + ^л
Че1
"02 " е2 ( c \
+ ^4 J gcos0dе +(a1 -с11 )/п + ^5 J gsin2edе + | b2 + -m° I = min
_ei _ Л ^ 4 _
Для классической обратной задачи выражения для множителей Лагранжа удалось получить в замкнутом виде [9]. Для обобщенной обратной задачи получим следующую систему линейных уравнений 5-го порядка:
Г N 0
0
12
s6 bs9
0 as4
s bs
^as3 bs9 as4 bs8 s11 у
Xo
Xa
VX5 у
=z
Здесь введены следующие обозначения:
1
c-г
-cA
V-c5 У
N
=z-
N
51 Z-, 2 2:
n=1 n - p
=z
n=1
(-2 - p2)
N
1
N
N
=z
=z-
1 n2 - 4pr 4 £ (n2 - p2 )(n2 - 4p2У 5 £ („2 - 4p2 )
-; s6 = N + p2s1; s7 = s1 + p2s2;
N 1
2 2 2 v ч 1 2 2
s8 = S3 + 4 p S4; S9 = N + 4 p S3; S10 = S3 + 4 p s5; su =^-y; S12 = b sw + as
n=1 n
p =
2n
S2-ег
d0 =- —; a = -4 p2 (cos ( 2S2 )-cos ( 2S1)); b = 2 p ( sin ( 2S2 )-sin ( 2S1));
Сл =
v1S12 + u1c12 .
s,n - na
Ui —
*0 12
1.
U1S12 v1c12 , (S122 + c122 ) p
v1 =-doC12 + ; c12=cosS2 -cosS1; s12 = sinS2 -sinS1;
С = -
12 12
e12 d0 + nb2 .
12
= cos (2S2 )- cos (2S1). Коэффициенты Фурье в выражении (10) находим
по формулам:
X,
Х-5
X5a
2 2 2 Л 2
n - p n - 4 p
X.
X 4 n
X 5n
2 2 2 2 „ 2
n - p n - 4 p
Решая систему уравнений, находим модифицирующую функцию £, а значит и функции Нт и ут. Снова по формулам (4) получаем решение обобщенной задачи для несжимаемой жидкости. При этом обеспечивается целевое распределение скорости по верхней поверхности профиля.
Для нелинейной задачи, когда скорость и (0) в расчетной плоскости в уравнении (2) не известна, как и в [9], будем решать обобщенную обратную задачу методом последовательных приближений.
Решение находим в ходе итерационного процесса. Считая функцию 5 (0) известной на к-й итерации, решаем прямую задачу для уравнений полного потенциала и находим скорость
(0). Уравнение (2) примет вид
sk+1 _ uk
(е)
d е u (s )
Интегрируя, находим промежуточное следующее приближение 5к+1, решаем систему уравнений (11) и определяем gk. Находим следующее приближение, удовлетворяющее условиям (8):
ds
k+1 _ UAk+1
ехр gk. И так до достижения критерия сходимости
k+1
< б, где б — малый
параметр. Во время итераций при решении прямой задачи задается не величина угла атаки, а значение параметра у = Г/Гт;п , а угол атаки определяется в процессе решения из условия Жуковского.
s
1
1
1
s
2
1
s
n
u
2. Рассмотрим приложение обобщенной обратной задачи к проектированию вертолетных профилей. При аэродинамическом проектировании вертолетных профилей обычно накладываются различные ограничения на величину ет0 профилей при использовании в различных компоновках вертолетных лопастей. Для определенности предположим, что требуется увеличить ет0 в сторону кабрирования. Например, это бывает необходимо после применения метода оптимизации передних кромок профилей [11]. Такое увеличение ет0 для вертолетных профилей обычно достигается за счет отгибания хвостовой пластины вверх. Это в свою очередь приводит к снижению величин су тах и Мкр. В нашем случае решается обобщенная обратная задача с сохранением
распределения давления по верхней поверхности. Поскольку в первом приближении сохраняется и целевой су (точно для Мот = 0), следует ожидать, что сутах модифицированного и исходного
профилей будут близки.
В качестве исходного профиля возьмем известный вертолетный профиль КЛСЛ-23012. В хвостовой части профиля вдоль его хорды установлена пластина. В качестве целевых возьмем распределения давлений, рассчитанные прямым методом [10] на трех характерных режимах обтекания: Мот = 0.4, су = 1.4 (режим сутах); Мот = 0.6, су = 0.6 (режим Ктах); Мот = 0.75,
су = 0 (режим Мкр).
Рис. 1. Распределения давлений и контура профилей (М = 0.4, су = 1.4)
-В-Є-Е*-X-BOOK-КЮ
Рис. 3. Распределения давлений (M = 0.6, су = 0.6)
Будем решать обобщенные обратные задачи при условии увеличения ст0 исходного профиля на 0.01. Обозначим полученные в результате решения профили следующим образом: N230^, N230^, N230^5, исходный профиль — N230, исходный профиль с отклоненной на 2° вверх хвостовой пластиной — N230
На рис. 1 приведены распределения ср и контуры профилей N230 и его модификации Nm4, полученные при решении обобщенной обратной задачи по данным прямой задачи при = 0.4, су = 1.4. На рис. 2—4 сравниваются распределения ср на трех исследуемых режимах для трех
модификаций исходного профиля и профиля с отклоненной пластиной. Как и следовало ожидать, семейство профилей, получаемое из решения трех обобщенных обратных задач на различных режимах обтекания, характеризуется распределением давления по верхней поверхности, близким к исходному профилю, а наибольшее отличие распределений давления по нижней поверхности наблюдается при = 0.75. На рис. 5 приведены распределения ср и контура профилей N230 и
его модификаций Nm4a, Nm4b и Nm4c, которые отличаются протяженностью зоны модифика-
ции, выраженной в долях суммарной длины дуги нижнего контура: 0.94, 0.97 и 0.99. Из рисунков видно, что полученные при решении обобщенной обратной задачи профили могут заметно отличаться формой нижнего контура и, соответственно, распределением давления по нему. Это позволяет на основе выбора режима обтекания (Mх, су ) и/или длины зоны модификации находить
оптимальное для заданной целевой функции решение обобщенной обратной задачи.
Коэффициенты продольного момента ст0 и максимальной подъемной силы суmax профилей N230, N230^, N230^ при Мот = 0.4 были рассчитаны по CFD модели ЯЛ№ [12]. По сравнению с исходным профилем были получены следующие величины приращений: Лст0 = 0.0119, Асутах =-0.0048 для профиля Ш30т4 и Лст0 = 0.0088, Асутах =-0.021 для
профиля N230 Расчетные данные подтверждают заявленное выше преимущество профиля
N230,
т4’
полученного в результате решения обобщенной обратной задачи, над профилем
N230^, полученным в результате отклонения хвостовой пластины вверх.
0Л -|V
Рис. 5. Распределения давлений и контура профилей (M = 0.4, cy = 1.4)
Таким образом, на основе решения обобщенной обратной задачи предложен метод модификации нижней поверхности профиля, позволяющий увеличить продольный момент профиля в сторону кабрирования без заметного снижения его максимальной несущей способности. Это позволяет использовать метод при аэродинамическом проектировании профилей, в частности для вертолетных лопастей.
ЛИТЕРАТУРА
1. Lighthill M. J. A new method of two-dimensional aerodynamic design // R&M 2112,
April 1945, ARC, England.
2. A r l i n g e r G. An exact method of two-dimensional airfoil design // TN 67, Oct., 1970,
Saab, Sweden.
3. Van Ingen J. L. A program for airfoil section design utilizing computer graphics //
AGARD short course notes, 1969.
4. Strand T. Designing airfoils in incompressible flow // J. Aircraft. 1973. V. 10, N 11.
5. Woods L. C. Airfoil design in two-dimensional subsonic compressible flow // R&M 2845, ARC, England, 1952.
6. Шагаев А. А. Построение контура профиля по заданному распределению давлений в сжимаемом потоке газа // Труды ЦАГИ. 1978, вып. 1925.
7. V o l p e G., Melnik R. E. The design of transonic airfoils by a well-posed inverse method // Intern. J. for numerical methods in engineering. 1986. Vol. 22.
8. Никольский А. А. Обратная задача для профиля при его потенциальном обтекании // Труды XIV Научно-технической конференции по аэродинамике и динамике винтокрылых летательных аппаратов, 1987.
9. Nikolsky A. A. Some aspects of helicopter airfoil design // Twenty first european rotorcraft forum. 1995. V. 2, N 17.
10. Bauer F., Garabedian P., Korn D. Supercritical wing sections III // Lecture Notes in Economics and Math. Syst. 1977. N 150.
11. Никольский А. А. Оптимизация передних кромок вертолетных профилей // Ученые записки ЦАГИ. Т. XXXIX, № 4, 2008.
12. Morrison J. H. A compressible Navier — Stokes solver with two-equation and Reynolds stress turbulence closure models // NASA CR-4440, May 1992.
Рукопись поступила 17/12012 г.