CC BY
74
УДК 629.735.33.014.16:533.662.6
ОПРЕДЕЛЕНИЕ АЭРОДИНАМИЧЕСКИХ ХАРАКТЕРИСТИК ПРОФИЛЯ ПО ПОЛЯРЕ КРЫЛА КОНЕЧНОГО РАЗМАХА
Б.Л. АРТАМОНОВ
Статья представлена доктором технических наук, профессором Ципенко В.Г.
Изложен численный метод расчета аэродинамических характеристик профиля на основе полученных из эксперимента аэродинамических характеристик прямоугольного крыла конечного удлинения. Задача решается методом последовательных приближений с использованием матрицы частных производных коэффициентов подъемной силы и аэродинамического сопротивления крыла по соответствующим характеристикам профиля.
Ключевые слова: численный метод, расчет аэродинамических характеристик, профиль.
В практике аэродинамических расчетов обычно необходимо иметь аэродинамические характеристики профиля на режимах досрывного и засрывного обтекания при различных числах Маха или Рейнольдса. Такие характеристики можно получить либо по результатам продувок дренированной модели профиля в плоской аэродинамической трубе, либо пересчетом аэродинамических характеристик модели крыла на бесконечное удлинение. Наиболее часто используется второй путь.
Пересчет аэродинамических характеристик крыла на бесконечное удлинение обычно базируется на методе Прандтля [1], который дает хорошие результаты в линейной области зависимости су(а). При расчете аэродинамических характеристик винтов этого диапазона недостаточно и необходимо иметь результаты в закритической области, причем как в зоне положительных, так и отрицательных углов атаки. Для решения этой задачи воспользуемся методом расчета аэродинамических характеристик крыла, основанным на его вихревой модели с непрерывно распределенной циркуляцией вихревого слоя по его размаху [2]. Решим обратную задачу: имея зависимости сук(фк) и Схк(фк), получим зависимости сур(ар) и схр(ар) в диапазоне углов ар, необходимом для расчета аэродинамических характеристик крыла.
Она осложняется тем, что при решении прямой задачи коэффициент подъемной силы и коэффициент сопротивления крыла зависят одновременно от сур(ар) и схр(ар), которые, как показано в [2], находятся под знаком интеграла по размаху крыла
1
°ук = ^ I \ур И С08 ^О) + схр (а) зЬ 1¥\г)Ъ(г)аг-
(о
сл= = \ I О) соз - сур (а) зЬ Дг) ж\г)Ь(г)с1г,
где '(7) = у]у2 - результирующий вектор скорости в сечении крыла; уу - вертикальный
компонент индуктивной скорости в сечении крыла; Р^) = Лгс1§[уу(2)/У] - угол притекания потока в сечении крыла; а^) = ф^) - Р^) - угол атаки результирующего вектора скорости в сечении крыла.
С математической точки зрения задача сводится к численному решению системы интегральных уравнений (1) относительно неизвестных функций сур(ар) и схр(ар) при заданных начальных условиях. Применим для ее решения метод последовательных приближений с использованием частных производных.
Суть предлагаемого метода заключается в том, что зависимости сур(ар) и схр(ар) последовательно выстраиваются относительно начальной точки ао, в которой коэффициенты подъемной силы и сопротивления профиля однозначно определяются характеристиками крыла.
Для плоского прямоугольного крыла имеют место следующие соотношения:
Cyp(ao) = Cyk(ao) = 0; Cxp(ao) = Cxk(ao). (2)
Предположим, что характеристики крыла заданы в некоторых дискретных точках сук = сук(фы), е^ = с^фи), 0 < i < М. Тогда характеристики профиля могут быть получены в виде дискретного набора точек сур(ар0, ^(а^), количество которых будет равно М.
Допустим, что зависимости сур(ар), с^а^ в диапазоне а0 < ар, < аpn уже найдены и аппроксимируются линейными сплайнами, проходящими через точки сур(ар0, е^а^), 0 < i < п. Требуется расширить область определения искомых зависимостей и определить характеристики профиля в следующей точке аp п+ь Cyp n+l(аp п+1) и Cxp n+l(аp п+1), исходя из условий Cyk п+1 = Cyk(фk п+1),
Cxk П+1 = Cxk(а фk П+1).
Будем искать характеристики профиля для угла атаки ар, п+1 = фk п+1. Поскольку на докрити-ческих режимах обтекания углы атаки в сечениях крыла не превосходят его угла установки по отношению к невозмущенному потоку, задание характеристик профиля линейным сплайном в диапазоне а0 < аp < аk п+1 является достаточным для вычисления по ним характеристик крыла при угле установки фk п+1.
Для нахождения сур п+^ар п+1) и cxp п+^ар п+1) воспользуемся методом последовательных приближений. Допустим, что после у-го приближения характеристики профиля на угле атаки ар п+1
равны с
урп+1
и с
(V)
крыла на фк„+1 равны схук'п+1 и с^к/п_1
хрп+1> а вычисленные по ориентированным на них сплайнам характеристики . Тогда невязка на этой итерации составит:
(3)
Найдем, на сколько должны измениться аэродинамические характеристики профиля сур п+1 и с^ п+1, чтобы невязка обратилась в ноль. Для этого вычислим по формулам (1) частные производные функций сук и Схк; по искомым переменным сур п+1 и cxp п+1 при фk = фk п:
д>) _>> - с • укп+1 — укп+1 укп+1;
лС ^ = с^ - с
хкп+1 хкп+1 хкп+1'
дс
ук
дс
ур п+1
дс
ук
дс.
<Рк=<Ркп
хр П+1
дс
хк
дс
<Рк=<Ркп
ур П+1
дс
хк
дс
<Рк=<Рк„
хр п+1
(4)
'Рк=<Ркп
и решим систему двух линейных уравнений
Дс
дс
ук
урп+1
с
Фк=Фкп
с
хрп+1 5схк
с
хрп+1
= дс<
укп 1
к кп
урп 1
с
к кп
хрп+1
с
хрп+1
= ДсМ
хкп 1
Фк=Фкп
относительно искомых поправок дсур п+1 и Дcxp п+1. Получим:
дс
хк
дс
Дс.
хр п+1
Дс,
(V)
укп+1
дс,
ук
дс
<Рк=<Ркп
хр П+1
Дс
хк п+1
<Рк=<Ркп
урп+1
Дс,
хр п+1
5сук дсЛ дСук 5схк
5сурп+! (Рк=(Ркп 5Схрп+1 <Рк=<Ркп 9 Схрп+1 4>к~- -~(Ркп 5Сурп+1 (Рк=(Ркп
&е{ v) хк п+1 <Рк=<Ркп ЗСхк &е{ v) укп+1
^урп+1 5Сурп+1 (Рк =<Ркп
^к
3Сурп+1 <Рк=<Ркп 3Схрп+1 (Рк=(Ркп д Схрп+1 =<Ркп 5сурп+1 <Рк=<Ркп
V
V
Следующее приближение искомых переменных формируем, используя общий для обеих переменных коэффициент релаксации ц:
Сурп+1 — Сурп+1 + п+1' Схрп+1 — С!рп+1 ^ ■ (6)
Отметим, что для обеспечения устойчивой апериодической сходимости итерационного процесса на всех углах атаки крыла, включая закритические, коэффициент релаксации должен лежать в диапазоне 0,4 < ц < 0,1. Приближения повторяются до тех пор, пока не удовлетворится условие окончания итерационного процесса:
Ас
W
у к п+1
Ас;
хк п+1
(7)
где с v, с х - абсолютные точности вычисления аэродинамических характеристик крыла, принятые равными sy = 0,01; sx = 0,001.
При удовлетворении этих условий переходят к расчету следующей точки поляры. По окончании расчета одной ветви кривых (при cy > 0) переходят к другой ветви (при cy < 0), начиная с той же исходной точки (2).
Работоспособность алгоритма, определяемого формулами (3) - (7), подтверждена расчетами аэродинамических характеристик профилей NACA 00-12 и NACA 230-10, испытания отсеков которых с X = 5 были проведены в ЦАГИ при различных числах Маха (0,3 < М < 0,9) в достаточно широком диапазоне углов установки (-5° < фк < 22°).
На рис. 1 показаны исходные аэродинамические характеристики отсека крыла cyk = cyk(9k,M), cxk = cxk(9k,M) с профилем NACA 00-12, включающие в себя как досрывную, так и засрывную области обтекания. Видно, что при больших числах Маха закон изменения коэффициента подъемной силы по углу атаки крыла существенно отличается от линейного, а максимальное значение коэффициента подъемной силы не превосходит 0,95 при М = 0,3. Поэтому применение методики пересчета, основанной на гипотезе линейности зависимости cy(a), было бы неправомерным.
^(ф^М) ^(ф^М)
Крыло
-М=0.3
- М=0.4 М=0.5 М=0.6
- М=0.7
- М=0.8
- М=0.85
- М=0.9
Крыло
-М=0.3 М=0.4 М=0.5 М=0.6 М=0.7 М=0.8 М=0.85 -М=0.9
-10
0
5
10
5
20
25
Рис. 1. Аэродинамические характеристики отсека крыла X = 5 с профилем NACA 00-12
На рис. 2 показаны аэродинамические характеристики профиля NACA 00-12 cyp = ^(а^М), cxp = ^(а^М), полученные расчетом по вышеописанному алгоритму. Анализ показывает, что они качественно повторяют характер исходных кривых, но имеют более высокие максимальные значения коэффициента подъемной силы cy max = 1,05 при М = 0,3, что соответствует данным для этого профиля, приведенным в [3].
На рис. 3 те же аэродинамические характеристики представлены в виде семейств поляр крыла с^^М) и профиля с^^М). Видно, что при всех числах Маха поляры профиля получены и в засрывной области обтекания, а точка ветвления поляр профиля смещена от аналогичной точки поляр крыла примерно на 0,1 в сторону меньших значений коэффициента сопротивления, что обусловлено отсутствием индуктивного сопротивления.
.ур(ар,М)
-хр(ОрМ)
10 - i 1 01 52 02
/-°2 ^0 4
—-0.6 -
Профиль
♦ М=0.3
• М=0.4
М=0.5
М=0.6
М=0.7
-•- М=0.8
—.— М=0.85
- М=0.9
Профиль
-М=0.3 М=0.4 М=0.5 М=0.6 М=0.7 М=0.8 М=0.85 -М=0.9
25
Рис. 2. Аэродинамические характеристики профиля NACA 00-12
Сук(СхкМ)
1.2 1
0.8 0.6 0 . 4 0.2 0 --0.2 ., i
-0 4 -Лг -0.6 -
Сур(СхрМ)
Профиль
-М=0.3 -М=0.4 М=0.5 М=0.6 -М=0.7 -М=0.8 - М=0.85 -М=0.9
Рис. 3. Поляры профиля и крыла с профилем NACA 00-12
Более наглядное отличие поляр профиля NACA 0012 от поляр крыла с тем же профилем следует из рис. 4, где они совмещены для одних и тех же чисел Маха. Там же треугольными маркерами нанесены поляры профиля, полученные на МВЗ им. М.Л. Миля пересчетом из тех же поляр крыла. Видно, что в досрывной области для всех чисел Маха они практически идентичны предлагаемым нами, а в засрывной области дают более высокие значения коэффициента подъемной силы. Причина этого обстоятельства заключается в том, что на МВЗ им. М.Л. Миля в досрывной и засрывной областях обтекания применялись разные методики пересчета, поскольку алгоритм, использовавшийся в досрывной области, за срывом не обеспечивал сходимости.
Предлагаемый в настоящей статье метод построен на единой схеме итерационного процесса в обеих областях, что и обусловило полученный результат. Расчеты показали, что если по характеристикам профиля, полученным на МВЗ им. М.Л. Миля, обратно рассчитать характеристики крыла, то в досрывной области они совпадут с исходными, а в засрывной области - дадут более высокие значения коэффициента подъемной силы, чем исходные. Полученные нами характеристики профиля этих недостатков не имеют.
Следует отметить, что при М = 0,3 - 0,4 в досрывной области при равных углах установки коэффициент подъемной силы профиля всегда больше, чем коэффициент подъемной силы крыла. Этот факт свидетельствует о том, что на крыле имеет место одновершинный закон циркуляции r(z), создающий индуктивное поле, уменьшающее углы атаки в сечениях крыла по сравнению с углом установки. На закритических углах установки крыла и закритических углах атаки профиля в сечениях это обстоятельство может привести не к уменьшению, а к увеличению коэффициента подъемной силы, поскольку зависимость cyp = сур(ар) имеет отрицательный градиент. Поэтому коэффициент подъемной силы крыла может оказаться больше коэффициента подъемной силы профиля, особенно в тех случаях, когда зависимость сур = сур(ар) в зоне срыва имеет резкое падение, что характерно для малых чисел Маха.
1.5
0.5
-0.5
1.5
0.5
-0.5
1.5
0.5
-0.5
ll 0 1 0.
М=0.3
03
- Крыло
Профиль
Пересчет МВЗ И-
04
1.5
0.5
-0.5
-1 1.5
0.5
-0.5
-1 1.5
0.5
-0.5
Крыло Профиль Пересчет МВЗ
М=0.6
03
- Крыло Профиль Пересчет МВЗ
~Г
Крыло Профиль Пересчет МВЗ
1.5
0.5
-0.5
М=0.85
0 W 0 1 0. 2 0. 3 0
• Крыло Профиль —ь— Пересчет МВЗ i
1.5
М=0.9
/
Ч 1 0 2 0 3 0.
Ф Крыло Профиль 1
0
0
4
2
0
0
4
0
0
0
2
0
4
0
0
4
4
0.5
0
0
4
4
0.5
Рис. 4. Сравнение поляр профиля NACA 00-12 и крыла с тем же профилем
Рассмотрим второй пример. На рис. 5 показаны исходные аэродинамические характеристики отсека крыла cyk = cyk(9k,M), cxk = cxk(9k,M) с профилем NACA 230-10, также включающие в себя как досрывную, так и засрывную области обтекания. В отличие от профиля NACA 0012,
закон изменения коэффициента подъемной силы по углу установки крыла при М = 0,3 и М = 0,4 имеет волнообразный характер, что обуславливает немонотонный характер зависимости Г^) и затрудняет сходимость итерационного процесса.
Сук(фк,М)
Cy
Схк(фк,М)
1.2 л
10 0 2 ( I 1 02 0 3
-0.6 -
Крыло
М=0.3
• М=0.4
М=0.5
-х- М=0.6
М=0.7
• М=0.8
—.— М=0.85
- М=0.9
0.45 У
&
//}
/¿7/ /
//У к
/
'Л ш(
"tro1 ■у*
Крыло
-М=0.3 -М=0.4 М=0.5 М=0.6 -М=0.7 -М=0.8 - М=0.85 -М=0.9
Cx
-10
0
10
20
30
Рис. 5. Аэродинамические характеристики крыла X = 5 с профилем NACA 230-10
Это обстоятельство потребовало некоторой корректировки описанного выше метода извлечения характеристик профиля из характеристик крыла путем формализации алгоритма вычисления коэффициента релаксации ц в формулах (6) в зависимости от знака и абсолютной величины частных производных (4). Это позволило обсчитывать и нелинейные зависимости
Сук(фк,М), Схк(фк,М).
На рис. 6 показаны аэродинамические характеристики профиля NACA 230-10 cyp = cyp(ap,M), cxp = cxp(ap,M), полученные расчетом по модернизированному алгоритму. Анализ показывает, что зависимости cyp(ap,M) качественно повторяют волнообразный характер исходных кривых сук(фк,М), но имеют более высокие максимальные значения коэффициента подъемной силы профиля сур max = 1,2 при М = 0,3 по отношению к коэффициенту подъемной силы крыла сук max = 1,04. Отметим, что полученные значения сур max соответствуют аэродинамическим характеристикам профиля NACA 230-10, приведенным в [3].
Cyp(ap,M)
Cxp(ap,M)
1.4 -|
г
J
0 /JL 0 1 1 2 0 3
-0.6 -
Профиль
♦ М=0.3
• М=0.4
М=0.5
М=0.6
-•- М=0.7
■ М=0.8
М=0.85
- М=0.9
-0т5-п
X
/ А I
Ж
jyy
/¿V /,
/ г
Профиль
-М=0.3 М=0.4 М=0.5 М=0.6 М=0.7
- М=0.8
- М=0.85
- М=0.9
Cx
-1
-10
0
10
20
30
Рис. 6. Аэродинамические характеристики профиля NACA 230-10
На рис. 7 поляры профиля и крыла совмещены для одних и тех же чисел Маха. Треугольными маркерами нанесены поляры профиля, полученные в ЦАГИ пересчетом из тех же поляр крыла. Видно, что для всех чисел Маха и в досрывной, и в засрывной областях зависимости практически идентичны предлагаемым нами, но приведены в диапазоне углов атаки не превосходящем 15°. Этот факт может быть связан как с отсутствием данных в рассматриваемой области, так и с недостаточной корректностью использования единой вихревой модели для расчета характеристик крыла на больших закритических углах атаки, где характер обтекания и структура вихревого следа имеют принципиальные отличия от докритических режимов.
Рис. 7. Сравнение поляр профиля и крыла с профилем NACA 230-10
Аналогичное сравнение было выполнено для зависимостей коэффициента подъемной силы от угла атаки. И здесь результаты расчета характеристик профиля хорошо совпали с данными ЦАГИ.
Отметим, что по данным ЦАГИ при угле установки 9k = 15° и М = 0,3 коэффициент подъемной силы крыла становится равным коэффициенту подъемной силы профиля. Это подтверждает сделанный выше вывод о том, что на засрывных режимах обтекания на крыле имеются участки, где угол атаки в сечении больше, чем угол установки.
Выводы
1. Математическая модель крыла конечного размаха, построенная на основе вихревой пелены с непрерывно распределенной по ее участкам погонной циркуляцией вихревого слоя, дает возможность получать расчетным путем аэродинамические характеристики профиля, основываясь на аэродинамических характеристиках крыла, взятых из эксперимента.
2. Метод последовательных приближений с использованием частных производных позволяет получать аэродинамические характеристики профиля из аэродинамических характеристик крыла конечного удлинения в области докритических и закритических режимов обтекания как при положительных, так и при отрицательных углах атаки.
3. Работоспособность алгоритма продемонстрирована в диапазоне углов атаки профиля -5° < а <25° и чисел Маха 0,3 < М < 0,9, а достоверность подтверждена сопоставлением полученных результатов с данными ЦАГИ и МВЗ им. М.Л. Миля.
ЛИТЕРАТУРА
1. Юрьев Б.Н. Крыло с изменяющейся по размаху циркуляцией: избранные труды. - М.: Изд-во АН СССР, 1961. - Т. 2.
2. Артамонов Б.Л. Вихревая модель крыла с непрерывно распределенной циркуляцией вихревого слоя // Научный Вестник МГТУ ГА. - 2014. - № 200. (В настоящем вестнике).
3. Миль М.Л. и др. Вертолеты // Расчет и проектирование // Аэродинамика. - М.: Машиностроение, 1966. - Т. 1.
DEFINITION OF AERODYNAMIC CHARACTERISTICS OF THE PROFILE UNDER THE DRAG POLAR OF THE FINITE WING
Artamonov B.L.
The numerical method of calculation of aerodynamic characteristics of a profile is stated on the basis of the experiment of aerodynamic characteristics of a rectangular wing of a final lengthening. The problem is solved by a method consecutive approximation with the use of a matrix of private derivative factors of elevating force and aerodynamic resistance of a wing under the appropriate characteristics of a profile.
Key words: numerical method, calculation of aerodynamic characteristics, profile.
Сведения об авторе
Артамонов Борис Лейзерович, 1947 г.р., окончил МАИ (1972), старший научный сотрудник, кандидат технических наук, заместитель заведующего кафедрой проектирования вертолетов МАИ (НИУ), автор более 180 научных работ, область научных интересов - аэромеханика винтокрылых летательных аппаратов вертикального взлета и посадки.
