______УЧЕНЫЕ ЗАПИСКИ Ц А Г И
Том XIX 19 8 8
№ 5
УДК 533.6.011.35:533.695
ПОСТРОЕНИЕ ПОВЕРХНОСТИ КРЫЛА В НЕВЯЗКОМ ОКОЛОЗВУКОВОМ ПОТОКЕ ПО ЗАДАННОМУ РАСПРЕДЕЛЕНИЮ ДАВЛЕНИЯ
А. Е. Осовский, В. И. Савицкий
Рассматривается задача построения формы поверхности крыла конечного размаха, имеющего заданное распределение давления и заданное распределение толщины задней кромки с учетом влияния фюзеляжа. Разработан численный алгоритм ее решения в рамках околозвуковой теории' малых возмущений. Приводятся примеры расчетов.
В последние годы численные методы невязкого околозвукового течения получили значительное развитие, в результате чего стало возможным их применение к проблеме проектирования формы аэродинамических тел с заданным распределением давления в плоском [1, 2] и трехмерном [3—5] потоке газа.
В работе [3] для решения обратной задачи аэродинамики крыла конечного размаха при околозвуковой скорости набегающего потока используется метод интегральных соотношений, а в работах [4, 5] развивается примененный в [2] подход, состоящий в замене граничных условий Неймана условиями Дирихле.
В настоящей статье алгоритм решения обратной задачи аэродинамики строится с использованием работы [6], где разработана численная программа и получено решение задачи обтекания комбинации крыла с фюзеляжем околозвуковым потоком газа. Для построения искомой поверхности крыла используется замена граничных условий Неймана условиями Дирихле, при этом для нахождения значений потенциала в плоскости крыла по заданному распределению давления используется численная процедура. Достижение заданного распределения толщины задней кромки осуществляется итерационно, внутри каждого цикла применяется метод релаксаций.
Уравнение изоэнтропийного течения невязкого газа записывается в рамках околозвуковой теории малых возмущений [6]:
-'4'Р хх + В (? у у + ^гг) — срг <р Х2 — 0, (1)
где
А = 1 — - (1 + 7) <р^1 + _1_ <р^ ;
В — 1 + (1 — ч) Мот ’
х, у, г — декартовы координаты (х — в направлении набегающего потока, г — по размаху крыла), — число Маха невозмущенного потока, у — показатель адиабаты, «р (х, у, г) — потенциал возмущения.
На поверхности крыла выполняется линеаризованное условие непротекания
*/|у=±о = -^т—“• (2)
здесь а — угол атаки крыла, функции г/+(*, г) и у~(х, г) задают координаты соответственно верхней и нижней поверхности крыла. Помимо этого потенциал удовлетворяет условиям непротекания на поверхности фюзеляжа, условиям затухания возмущений на бесконечности и условиям на поверхности вихревой пелены, сходящей с задней кромки
крыла [6]. Коэффициент давления, на поверхности крыла вычисляется по формуле, учитывающей члены, квадратичные относительно производных потенциала возмущений
еР 1у = ±0 = - 2?х - (I - ЛО - ср2 -<р2. (3)
Для удобства проведения расчетов осуществляется переход к новой неортогональной системе координат 0£г||, позволяющей отобразить бесконечное физическое пространство на конечную расчетную область К I < 1, 1 > ^ > 0, 1 с | < I. Внутри этой области вводится равномерная расчетная сетка и для вычисления поля потенциала возмущения применяется метод верхней последовательной релаксации вдоль линий.
Рассмотрим задачу нахождения функций г/+(х, г) и у~{х, г), для которых выполняется соотношение
ср(х• У)1у = ±о = ср (*■ г>- (4)
где с+ (х, г) и с~(лг, г)—функции, задающие требуемое распределение давления со-ответственно по верхней и нижней поверхностям крыла.
Искомая форма поверхности крыла должна иметь заданную толщину задней кромки
/(*)=/*(*), |
Нг) = У+ К. в, г) — у~(х3 к, г), ) (5)
где х3. к (г)—координата задней кромки, 1*(г)—задаваемое распределение толщины задней кромки.
Данная задача вообще говоря, не является корректной, так как может не существовать формы поверхности крыла, для которой при заданном числе Мт выполняется соотношение (4). В работе [7] для плоского дозвукового потока получены условия разрешимости задачи, накладываемые на задаваемое распределение давления, а для рассматриваемого случая подобных условий нет. В связи с этим обычно формулируют обратную задачу аэродинамики, как задачу корректировки распределения давления на исходном крыле путем малой модификации формы его поверхности.
Будем искать распределение потенциала в плоскости крыла, используя формулы (3), (4). Замена этих соотношений их конечно-разностными аналогами приводит к системе квадратичных алгебраических уравнений, разрешая которые последовательно относительно потенциала возмущений получим
-а‘ЦУ (ар)2 - аЦк4к
Чцк =---------------уь-----------’ г=,гп. к> —> гз. к, } —/ф + Ь •"> ■/*> (6)
а\
здесь индекс г'п. к соответствует точкам на передней кромке крыла, 13. к — на задней кромке, /ф — задает граничные узлы, попадающие внутрь фюзеляжа, д — узлы на конце крыла, а индекс к соответствует точке на поверхности крыла. В выражения для входят значения коэффициентов преобразования координат в узлах сетки и известные значения потенциала возмущений в прилежащих узлах, причем при /=/ф + 1 для их определения дополнительно привлекается условие непротекания на поверхности фюзеляжа. Значения потенциала на передней кромке крыла считаются известными, они выбираются из условия заданной толщины задней кромки сечений крыла. В выражения для коэффициентов в формуле (6) входит также квадрат производной Фу в узлах сетки на поверхности крыла. Согласно околозвуковой теории малых возмущений член <Ру в (3) мал по сравнению с (р*, поэтому при расчете значений потенциала возмущений на поверхности крыла первоначально задается величина ц>у, рассчитанная для исходной геометрии крыла, а впоследствии она уточняется.
Найденное распределение потенциала возмущений на поверхности крыла позволяет решать задачу расчета околозвукового обтекания комбинации крыло — фюзеляж, заменяя условие непротекания (2) условием Дирихле
у(х, у=±0, г) = у±(х, г), (7)
где ф+(я, г) и ф~(х, г)—функции, значения которых вычисляются согласно (6). Для решения задачи, как и в предыдущем случае, используется метод верхней последовательной релаксации по линиям. Проведенные численные расчеты показали, что при осуществлении релаксационной процедуры граничные условия Дирихле дают более сильные возмущения, чем граничные условия Неймана. В связи с этим, во избежание появления численной неустойчивости необходимо специально подбирать оптимальную величину коэффициента верхней релаксации. В работе [6] коэффициент верхней релаксации £2 переменный по узлам сетки
., 1 . + ^ . 1 Л V х1+~2] х1~Т])
2 =-------*---------^-------> (8)
ЪхЦ ? !+-£-+!
х!-т1 “
ш
1,0 -
;о”
20°
ж
ЧП°
Рис. 1
где — коэффициент преобразования системы координат, а м — постоянное значение коэффициента релаксации на равномерной расчетной сетке в физическом пространстве.
В настоящей статье определялась оптимальная величина коэффициента ш. Как показали расчеты, она зависит главным образом от формы крыла в плане, и особенно от угла стреловидности. На рис. 1 приведен пример изменения оптимальной величины коэффициента релаксации при изменении угла стреловидности для комбинации крыла с фюзеляжем.
Поле значений потенциала возмущений, найденное в результате решения задачи Дирихле, используется для нахождения наклона поверхности крыла. Для этого уравнение (1) записывается в точке на.поверхности крыла, причем при вычислении производной используется граничное условие (2). Это приводит к выражению
ЗуМ _+ ёхкУцк±\ + Ш%кЧцк + .
)Ч ~ 2ВМ
дх
(9)
Л] к
здесь Ц'т определяются коэффициентами преобразования системы координат, взятыми в узлах сетки.
По вычисленному согласно (9) распределению наклона поверхности строится форма профилей в сечениях крыла. Для этого применяется кубическая сплайн-функция, которая на каждом интервале между двумя расчетными узлами задает координаты поверхности по формуле
у(х) =(-&.) (х - х.) + а,
дх и
(х-х& щ (х - х,у
6 Лг
+
+
+ °/+1
(х
■*/)3 + “>1 (X —
6Л;
-Мг. х1+1>х>х1г
(10)
где XI — узлы расчетной сетки, в точках известен наклон поверхности ду ' - х>
(—V
\ дх /г
Лг = х,
<+1
Для нахождения коэффи-
циентов и <1-1 решается система линейных уравнений с тридиагональной матрицей.
Задавая различные распределения потенциала на передней кромке крыла, получим при неизменном задаваемом распределении давления различные граничные условия в задаче Дирихле и, как следствие решения обратной задачи, дающие разные толщины задней кромки в сечениях крыла. Поэтому процедура построения заданного распределения толщины задней кромки сводится к нахождению потенциала на передней кромке крыла, что осуществляется итерационным способом. На каждой итерации поправка значений потенциала находится из уравнений
’к
2 АЧ А,?Ит = 1* (*) —1 (**)• АЧ = ’
1Чф+1
I = ^и. к> 2* — Уф 1, ..., Ук,
где индекс т соответствует точке на поверхности крыла. Коэффициенты матрицы Ац находятся численно по формуле
АЧ = ------>
(И)
(12)
где йу?—малое изменение уут, й/// — соответствующее ей приращение / (гг) после N шагов релаксационной процедуры, /* — приращение при отсутствии возмущений. Итерационный процесс (11) и (12) хорошо сходится уже при Л/=3.
Проверка работоспособности численного метода решения обратной задачи аэродинамики проводится для комбинации крыло — фюзеляж, типичной для современного пассажирского самолета. При решении обратной задачи существенным является выбор задаваемого распределения давления, что представляет собой самостоятельную сложную задачу. В первом приведенном примере расчета используется концепция, согласно которой линии постоянного давления должны быть прямыми и иметь угол наклона, равный углу стреловидности крыла. Решается задача построения формы поверхности крыла, имеющего одинаковую эпюру давления во всех сечениях по размаху. На рис. 2—4 приведены эпюры давления, расположение изобар и геометрические характеристики построенного крыла. Сравнение распределения давления на построенном крыле с заданным распределением давления, а также расположение линий постоянного давления показывают, что поставленная цель достигнута. По мере приближения к борту фюзеляжа крутка резко нарастает, что согласуется с результатами, полученными для дозвукового потока в рамках линейной теории [7].
Во втором примере рассматривается околозвуковое обтекание комбинации крыло— фюзеляж, типичной для современного пассажирского самолета, при Моо = 0,8 и су = 0,5. С помощью разработанного метода решается смешанная задача аэродинамики: распределение давления на верхней поверхности крыла от борта фюзеляжа до мес-
Ср на построенном крыле.--------------Ср заданное
Рис. 2
п
Верхняя поверхность крыла Нижняя поверхность крыла
Рис. 3
та излома задней кромки модифицировано так, чтобы устранить наличие двух областей сжатия в окрестности х=0,2 и л:—0,6; на этой части крыла определяется форма его поверхности. Модифицированное распределение давления выбирается таким образом, чтобы относительная толщина профилей в сечениях крыла сохранялась неизменной. Для этого, по возможности, поддерживается заданный исходным распределением давления уровень средней линии эпюры (Ср -)- с~)/2, который. согласно приближенной теории тонкого профиля определяет его толщину. Другая часть поверхности от места излома до конца крыла остается неизменной — на ней задается условие непротекания (2).
На рис. 5 показано сравнение распределения давления в сечениях построенного в результате расчета крыла с исходным и модифицированным (используемым в качестве входных данных для решения задачи) распределением давления.
На рис. 5 также представлена форма профилей исходного и построенного крыла. Сравнение приведенных результатов показывает, что желаемого эффекта удалось достигнуть за счет малых изменений формы профилей сечений крыла. На рис. 6 пока- ' заны изменения относительной толщины и крутки, полученные в результате применения предложенной методики решения смешанной задачи обтекания крыла. Распределение толщин профилей сечений крыла по его размаху практически осталось неизменным. Крутка построенного крыла в сечениях 0,1<г<0,4 меньше, чем крутка исходного крыла, при этом величина Ае не превышает 1°. Проверка полученных результатов с помощью программы расчета околозвукового обтекания крыла на основе численного интегрирования полного уравнения для потенциала [10] показала, что коэффициент сопротивления построенного крыла при постоянном значении коэффициента подъемной силы меньше на 8%, чем коэффициент сопротивления исходного крыла. Этот результат объясняется уменьшением волновой составляющей сопротивления в случае расчета обтекания крыла, построенного с помощью изложенного метода.
------исходное крыло;-----------ср модифицированное; --------построенное крыли
Рис. 5
ЛИТЕРАТУРА
1. Шагаев А. А. Построение контура профиля по заданному распределению давлений в сжимаемом потоке газа. — Труды ЦАГИ, 1978, вып. 1925.
2. Т г a n е n Т. L. A rapid computer aided transonic airfoil design method. — AIAA Paper 74-501, 1974.
3. Takanashi S. Iterative three-dimensional transonic wing using integral equations. — Journal of Aircraft, 1985, N 8.
4. H e n n e P. A. Inverse transonic wing design method. —
Journal of Aircraft, 1981, vol. 18, N 2.
5. S h a n k a г V., Malmuth N. D., Cole J. D. Computational
transonic inverse procedure for wing design. — AIAA J., 1982, vol. 20, N 8.
6. Савицкий В. И. Расчет невязкого околозвукового обтекания крыла с фюзеляжем. — Ученые записки ЦАГИ, 1983, т. 14, № 2.
7. Li gh thill М. J. A new method of two-dimensional aerodynamic design. — R&M 2112, April, 1945, ARC.
8. Рыжкова М. В., Серебрийский Я. М. Теоретическое исследование формы тонкого стреловидного крыла с постоянным распределением нагрузки по размаху. — Технические отчеты ЦАГИ, 1953.
9. Седов Л. И. Плоские задачи гидродинамики, аэродинамики.— М.: Наука, 1980.
10. Jameson A., Caughey D. A. Numerical calculation of the transonic flow past a swept wing. — N.-Y. University ERDA Report COO 3077-140, 1977.
Рукопись поступила 14/1V 1986 г.