Vol. 19, № 01, 2016
Civil Avition High TECHNOLOGIES
УДК 533.6.013.42
МЕТОД ДИСКРЕТНЫХ ОСОБЕННОСТЕЙ С МОДИФИЦИРОВАННЫМИ КВАДРАТУРНЫМИ ФОРМУЛАМИ В ЗАДАЧАХ АЭРОУПРУГОСТИ ЛЕТАТЕЛЬНЫХ АППАРАТОВ
В.В. ОВЧИННИКОВ, Ю.В. ПЕТРОВ, В.М. ПОПОВ
Рассматривается возможность применения модифицированных квадратурных формул метода дискретных вихрей при решении задач аэроупругости. Приводятся результаты расчетов для ряда тестовых задач, подтверждающие эффективность предлагаемых формул.
Ключевые слова: аэроупругость, дискретные особенности, квадратурные формулы, расчетные методы.
Развитие компьютерной техники позволило ставить и численно решать сложные задачи динамической аэроупругости с учетом нелинейности аэродинамических характеристик летательного аппарата. При формировании методов и построении методик расчета аэроупругой устойчивости и динамического нагружения летательных аппаратов с использованием линейных и нелинейных аэродинамических моделей необходимо особое внимание уделять вопросам корректного дискретного представления поверхности летательных аппаратов, выбору рациональных параметров разбиения, регуляризации решений систем линейных алгебраических уравнений (СЛАУ), к которым редуцируются соответствующие граничные интегральные уравнения. Данный вопрос очень важен для рассмотрения именно аэроупругих задач, так как при их решении первостепенное значение приобретает корректный расчет распределенных аэродинамических характеристик деформируемого тела, а дискретизация поверхности упругого летательного аппарата изменяется на каждом временном шаге численного расчета.
Настоящая статья посвящена анализу квадратурных формул для интегралов с сильной особенностью, которые встречаются при численном решении задач обтекания деформируемого тела методом дискретных вихрей (МДВ), и построению квадратур, слабо чувствительных к выбору алгоритмов создания вихревых схем и расстановки контрольных точек [1 - 3].
Известно, что задача обтекания упругого тела при малых дозвуковых скоростях сводится к задаче Неймана для уравнения Лапласа относительно потенциала скорости. Эта задача может быть приведена к сингулярному интегральному уравнению [1, 4] относительно интенсивности двойного слоя g с ядром н(г, |) вида:
J g(| М?, I )ds = ОД, (1)
S
которое выражает выполнение условия непротекания газа через поверхность тела в точке с радиус-вектором г (под f (г) здесь понимается нормальная составляющая скорости внешнего потока с учетом деформаций поверхности).
Сформулируем основные предположения и допущения относительно формы поверхности летательного аппарата (ЛА), свойств искомого решения и алгоритма схематизации поверхности.
П 1. Пусть S - ограниченная замкнутая поверхность. Положим, что орт внешней нормали n| существует и непрерывен во всех точках поверхности.
П 2. Диаметром поверхности будем называть величину DS = max г — I.
S (r,I)cS
Civil Avition High TECHNOLOGIES
Vol. 19, № 01, 2016
П 3. Пусть замкнутая поверхность Б является строго выпуклой, т.е. она может иметь не более двух точек пересечения с произвольной прямой.
П 4. Положим, что интегральное уравнение (1) имеет точное решение Б, удо-
влетворяющее условиям:
g(I)<Cl Vg(I)
< C
g(r) - g(I) - (V g(I),r -I)
< m.
r-I
2
где С1, С2, шё - некоторые положительные постоянные.
П 5. Фиксируем малое значение шага разбиения И и положим, что поверхность Б можно разбить на N непересекающихся областей, в каждой из которых выберем контрольную точку Я; так, чтобы выполнялись следующие 4 условия:
N • s=и sj, sj n sv=d sj n д sv при j , j=i • V? = г -> С0И УУ = 1,N, С0 - конечная константа,
• max r - I < h, Vv = 1, N, {r,I}<=Sv • ЗС3 > 0:Б = C3Nh2.
Для интегралов с сильной особенностью рассмотрим квадратурную формулу, учитывающую возможность линейного изменения интенсивности двойного слоя в пределах каждой из ячеек разбиения, т.е. для интеграла, понимаемого в смысле конечного значения по Адамару, вида:
1(8) = |)Е( г,
выберем
_ N ( ^ ^ _ _ ^ ^
квадратуру ЦГ) = £ g(Rv ) JS(r, I)ds + (Vg(Rv ), J(I- Rv )S(r, I)ds)
v=1
где
S(r, I) = (n, Vq>) =
(nr,n|)(r-I)2 -3(nr,r-|)(П|,r-I)
I'
?-I
и вектор Г = {g(Rv )} v = 1, N.
При выполнении этих условий справедлива Теорема 1. При выполнении условий П 1 -П 5 для произвольной контрольной точки замкнутой поверхности г = справедлива оценка
( П 2 ^ П
5 = |l(g)-Х(Г)| < 2nmg — + 4 h(l + O(h)).
C
V^ 0 у
Сделаем еще одно предположение.
П 6. Пусть поверхность такова, что в окрестности каждой из контрольных точек можно выбрать местную ортогональную систему координат ОХУ (с центром в контрольной точке) и составить матрицы численного дифференцирования Бх и Бу такие, что выполнены условия:
{V ^ (- яг
< mdh / 2,
{V ^ ( ^
Dy Г
=1,N y
< mdh / 2.
v
v
5
где ша > 0 - число, не зависящее от параметров дискретного представления поверхности.
Vol. 19, № 01, 2016
Civil Avition High TECHNOLOGIES
В контрольных точках г определим векторы:
WW1(r) = I J E(r, |)ds[ , Л1Х(У) = J J (I - ¿v )S(r, I)ds[ , Х(Г) = (\^ + A1XDX + A^Dy )г.
v=1,N
x(y) v=1,N
Тогда выполняется Теорема 2. При вы полнении условий П1 - П6 для произвольной контрольной точки замкнутой поверхности г = справедлива оценка
5 = |l(g) -ЦГ )|
<
2m,
g
( _2 ^ + 4
V Co
Л
+
mdn
у
2
2f 1 , 1 » — + ln—
C C
V^ o ^o
n h(1 + O(h)).
Аналогичные утверждения справедливы в рамках указанных выше предположений и для тонких поверхностей [1].
Поскольку ограничения на поведение искомого решения в рамках описанных выше исследований были достаточно жесткие, следует ожидать их нарушения на практике. Поэтому вторая часть исследований была посвящена оценке реальных преимуществ, которые предоставляет использование модифицированных квадратур. Естественно, для предварительного анализа были выбраны простейшие случаи, когда легко выявить причинно-следственные связи в наблюдаемых эффектах. Анализ проводился при численном решении как двумерных, так и трехмерных задач. Были отмечены следующие общие положительные черты (преимущества) предлагаемого модифицированного метода: во-первых, повышенная устойчивость результатов моделирования к высокой степени неравномерности разбиения поверхности (допустимо десятикратное скачкообразное изменение линейного масштаба панелей разбиения); во-вторых, повышенная устойчивость результатов моделирования к алгоритмам выбора положения контрольных точек (для равномерного разбиения квадратной пластины на квадратные панели допустима расстановка контрольных точек с использованием датчика случайных чисел внутри круга радиуса 0,45а с центром в геометрическом центре панелей, где а - сторона панели); в-третьих, модифицированный метод позволяет более устойчиво рассчитывать характеристики обтекания тонких тел (таких, у которых толщина сопоставима с шагом разбиения И), в-четвертых, использование модифицированных квадратур оправдано при решении не только линейных, но и нелинейных задач. Собственно "нелинейные" алгоритмы являются более ресурсоемкими по сравнению с "линейными", и исследователи, как правило, пытаются подобрать и использовать наиболее простую вихревую схему, "ухватывающую" основные детали компоновки. В частности, появился даже термин "настройка модели", означающий подбор схемы с минимальным числом особенностей, позволяющей, тем не менее, получать суммарные аэродинамические характеристики с заданной точностью (например, для дальнейшего совместного решения задач аэродинамики и динамики полета). Как показали исследования простейших компоновок (крылья простой формы в плане), модифицированные квадратуры позволяют на одном и том же разбиении с малым числом особенностей (К = 25 ■ 64) повысить точность определения суммарных аэродинамических коэффициентов на 5 % ■ 15 %.
Относительно свойств матрицы скосов (задача Неймана), определяемой интегралом
wjv = J
Sv
rA(?j, !v)
ds, A
R;
- Г
jv
: A(fj, Iv )
3(Пг,г-|)(П | ,r-I) - (Пг,П | )(r-I)
I-
установлено следующее.
v
v
2
2
3
Научный Вестник МГТУ ГА_Том 19, № 01, 2016
Civil Avition High TECHNOLOGIES Vol. 19, № 01, 2016
П 7. Пусть гладкая замкнутая строго выпуклая поверхность такова, что AjV > 2CA > 0, тогда такую поверхность будем называть правильной (здесь 2CA - некоторая положительная константа, зависящая от формы поверхности).
Теорема 3. Для замкнутых поверхностей при выполнении предположений относительно матрицы скосов справедливы следующие утверждения:
1. Для произвольной поверхности матрица имеет неполный ранг. Далее для правильных поверхностей:
2. Матрица W = (w;j) имеет ранг, равный N - 1, а ранг матрицы Wp системы, расширен-
N
ной условием ^Г = Cw, равен N. j=i
3. Норма обратной матрицы скосов Wu 1 для внешней задачи ограничена и не превосхо-
-i
дит величины о : о = min
V=1,N-1
wVN
, где Wu = {wV ■ }v=1,n-1 .
j=1,N-1
4. Все элементы обратной матрицы скосов 1 - положительные числа.
5. Регуляризированная матрица 'г имеет полный ранг, Г = Ь, где ^ = w ^, vN = 1, V= 1,К, j = 1,К -1.
6. Регуляризированная иным способом матрица скосов 'а имеет полный ранг, ее элементы вычисляются по формуле w = w ^ +а, V,j = 1,^ где а Ф 0 - произвольная регуляри-
зирующая константа.
: w
7. Матрица ' может быть приведена к матрице полного ранга wSv. = w V. +аБ.,
V,j = 1,N.
2
8. Если выбрать а = —Аз0, то
CAC
d3
Ws -1
3
<- S
D
SC Л ^ Г
Б
По сравнению с известными [2] алгоритмами построения невырожденных матриц СЛАУ МДВ для замкнутых поверхностей, предложенный в п. 7, обладает рядом преимуществ. С одной стороны, его использование приводит к наличию меньшей скорости протекания замкнутой поверхности тела и не требует выделения циркуляции одной из рамок в качестве регуляризи-рующей переменной (в соответствии с п. 5), с другой стороны, он не требует применения метода наименьших квадратов для определения значений циркуляций из решения расширенной СЛАУ (в соответствии с п. 2).
Если поверхность тела замкнутая, то регуляризированная СЛАУ при использовании модифицированных квадратур имеет вид:
НТ = b , S'j = {w + AxDx + AyDy + AzDz Jij + аSj, а *
В первую очередь, методическим исследованиям были подвергнуты алгоритмы, реализующие предложенные выше модификации квадратурных формул для интегралов с сильной особенностью. Изучалось влияние особенностей вихревых схем (неравномерность меры дискретности вдоль моделируемой поверхности, различные алгоритмы выбора контрольных точек) на точность результатов по суммарным и распределенным аэродинамическим характеристикам. Например, при решении линейной стационарной задачи для прямого крыла (А = 1) анализировалось влияние различных способов схематизации (рис. 1). На рис. 2-3 кривые (Series 3m - Series 5m) соответствуют случаю использования модифицированных квадратур, кривые (Series 3 - Series 5) - базовому МДВ, расчетная точка
Том 19, № 01, 2016_Научный Вестник МГТУ ГА
Vol. 19, № 01, 2016 Civil Avition High TECHNOLOGIES
(Series 2) соответствует данным [5]. Использование модифицированных квадратур приводит к существенному уточнению результатов моделирования в широком диапазоне изменения числа особенностей.
а - разбиение 1 типа (Series 1)
б - разбиение 3 типа (Series 3)
в - разбиение 4 типа (Series 4) г - разбиение 5 типа (Series 5)
Рис. 1. Анализируемые вихревые схемы (модельная задача)
Рис. 2. Изменение коэффициента подъемной силы по углу атаки Су с увеличением числа особенностей N
Civil Avition High TECHNOLOGIES
Vol. 19, № 01, 2016
m„
0.7
0.65
0.6
0.55
0.5
0.45
- Series 1
—1— Series 2
--Ж-- Series 3
■■в— Series 4
--Х— Series 5
Series 3m
-А— Series 4m
-н— Series 5 m
0 200 400 600 800 1000 N
Рис. 3. Изменение коэффициента момента тангажа по углу атаки m а с увеличением числа особенностей N
Нелинейные математические модели аэродинамики включают в себя алгоритмы регуляризации СЛАУ МДВ при наличии в схематизации поверхности ЛА замкнутых поверхностей. Рассмотрим следующий модельный пример. По результатам решения задачи о бесциркуляционном обтекании сферы единичного радиуса были вычислены безразмерные величины среднеквадратического отклонения давления от точного решения в контрольных точках для трех случаев (рис.4): решение строилось с использованием изложенного выше модифицированного алгоритма регуляризации (является базой для обезразмеривания), ре-гуляризирующей переменной, а также с использованием алгоритма на базе метода наименьших квадратов.
Предложенный алгоритм в несколько раз снижает погрешности вычисления распределенных характеристик, что позволяет более точно вычислять значение обобщенных аэродинамических ^ В """ """" N сил и, следовательно, адекватно анализировать ди- Рис. 4. Сравнение точности численного решения
5
6 5 4 3 2 1 О
M
/
г
Per, п.перем енная
Ы вдиф.м ?тод
намику упругих колебаний элементов летательного аппарата. В
в зависимости от числа особенностей N при использовании различных способов регуляризации (модельная задача)
качестве примера использования предложенного алгоритма регуляризации при решении нестационарной нелинейной задачи обтекания предлагаются результаты расчета (вихревой след за самолетом на рис. 5 и зависимости коэффициентов подъемной силы от безразмерного времени счета на рис. 6) компоновки самолета типа Як-130 на двух углах атаки а = 10° и а = 20°.
Vol. 19, № 01, 2016
Civil Avition High TECHNOLOGIES
Рис. 5. Пример расчета обтекания самолета типа Як-130 с использованием предложенного алгоритма регуляризации
Суа 1.4
1.2
0,8
0.6
0,4
0,2
0,0
A
V
^ 1 ■ ■ ■ ■
—■—œ=10 —•—a=20
5' X
Рис. 6. Изменение коэффициента подъемной силы по безразмерному времени при обтекании самолета типа Як-130
Таким образом, изложенные выше модификации МДВ с замкнутыми вихревыми рамками (модифицированные квадратурные формулы и предложенный алгоритм регуляризации) наряду с рациональными алгоритмами [3] вычисления скоростей в контрольных точках и узлах вихревой сетки позволяют существенно улучшить качество численных решений линейных и нелинейных задач аэроупругости современных ЛА при незначительном усложнении вычислительных алгоритмов.
ЛИТЕРАТУРА
1. [Попов В.М.| Математическое моделирование процессов взаимодействия летательного аппарата с внешними полями и разработка универсальных вычислительных процедур комплексного анализа аэродинамических компоновок. Дисс. на соискание уч. степ. докт. техн. наук. - М.: ВВИА им. проф. Н.Е. Жуковского. 2001. - 496 с.
2. Дворак A.B. Невырожденность матрицы метода дискретных вихрей в задачах пространственного обтекания. // Труды ВВИА им. проф. Н.Е.Жуковского. 1986. Вып. 1313. С. 441-453.
3. Ненашев A.C., Попов В.М. Использование формул приближенных вычислений при решении инженерных задач методом дискретных вихрей. // Научный вестник МЕТУ ЕА. 2008. № 125. С.102-105.
4. Морозов В.И., Овчинников В.В. Основы аэроупругости элементов боевых летательных аппаратов: монография. - М.: ВВИА им. проф. Н.Е. Жуковского. 2007. - 128 с.
5. Белоцерковский С.М., Моисеев Е.М., Табачников В.Е. Атлас нестационарных аэродинамических характеристик крыльев различной формы в плане. - М.: Изд. ЦАЕИ. 1959. - 478 с.
THE DISCRETE VORTEXES METHOD WITH THE IMPROVED QUADRATURE FORMULAE IN THE AEROELASTICITY PROBLEMS
Ovchinnicov V.V., Petrov Y.V., |Popov V.M.
The possibility of the improved quadrature formulae of the discrete vortexes method used in the aeroelasticity problems is considered. Some of the calculation results confirm the efficiency of these formulae are carried out.
Key words: aeroelasticity, discrete features, quadrature formulas, calculation methods.
Civil Avition High TECHNOLOGIES
Vol. 19, № 01, 2016
REFERENCES
1. |Popov V.M.| Matematicheskoe modelirovanie processov vzaimodeyjstviya letateljnogo ap-parata s vneshnimi polyami i razrabotka universaljnihkh vihchisliteljnihkh procedur kompleksnogo analiza aehrodinamicheskikh komponovok [Mathematical modeling of the interaction of the aircraft with external fields and the development of universal computational procedures comprehensive analysis of the aerodynamic configurations. Diss.] Diss. na soiskanie uch. step. dokt. tekhn. nauk. Moscow. VVIA im. prof. N.E. Zhukovskogo. 2001. 496 p. (In Russian)
2. Dvorak A.V. Nevihrozhdennostj matricih metoda diskretnihkh vikhreyj v zadachakh pros-transtvennogo obtekaniya [The nonsingularity of the matrix of the method of discrete vortices in problems of space flowing around. Proceedings of the Air Force Engineering Academy named after prof. Zhukovsky] Moscow. Trudih VVIA im. prof. N.E. Zhukovskogo. 1986. Iss. 1313. PP. 441 - 453. (In Russian)
3. Nenashev A.S., [Popov V.M.| Ispoljzovanie formul priblizhennihkh vihchisleniyj pri reshenii inzhenernihkh zadach metodom diskretnihkh vikhreyj [Using the formulas of approximate calculations in solving engineering problems the method of discrete vortices. Scientific Bulletin MSTUCA] Moscow. Nauchnihyj vestnik MGTU GA, seriya "Aehromekhanika i prochnostj". 2008. Iss.125. PP. 102 -105. (In Russian)
4. Morozov V.I., Ovchinnikov V.V. Osnovih aehrouprugosti ehlementov boevihkh letateljnihkh apparatov [Aeroelasticity basics of combat aircrafts elements]: monografiya. Moscow. VVIA im. prof. N.E. Zhukovskogo. 2007. 128 p. (In Russian)
5. Belocerkovskiyj S.M., Moiseev E.M., Tabachnikov V.G. Atlas nestacionarnihkh aehrodinamicheskikh kharak-teristik krihljev razlichnoyj formih v plane. [Atlas unsteady aerodynamic characteristics of wings of different shapes in the plan] Moscow. Izd. CAGI. 1959. 478 p. (In Russian)
СВЕДЕНИЯ ОБ АВТОРАХ
Овчинников Валерий Валерьевич, 1970 г.р., окончил МИФИ (1993), доктор технических наук, профессор, начальник кафедры механики и инженерной графики Академии ЕПС МЧС России, автор более 80 научных работ, область научных интересов - аэроупругость и прочность конструкций, электронный адрес: [email protected].
Петров Юрий Владимирович, 1956 г.р., окончил Рижское ВВАИУ (1978), доктор технических наук, профессор, заведующий кафедрой технической механики и инженерной графики МЕТУ ЕА, автор более 85 научных работ, область научных интересов - аэроупругость и прочность летательных аппаратов, электронный адрес: [email protected].
Попов Виктор Михайлович, 1962 г.р., окончил МФТИ (1985) и ВВИА им. Н Е. Жуковского, доктор технических наук, доцент, профессор МЕТУ им. Н.Э. Баумана, автор более 80 научных работ, область научных интересов - вычислительная аэрогидродинамика.