CC BY
45
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
То м XV
198 4
№ 3
УДК 629.735.33.015.4 : 533.6.013.422
КОМБИНИРОВАННЫЙ МЕТОД РАСЧЕТА АЭРОДИНАМИЧЕСКИХ СИЛ НА КОЛЕБЛЮЩЕМСЯ ЛЕТАТЕЛЬНОМ АППАРАТЕ В СВЕРХЗВУКОВОМ
ПОТОКЕ
Рассчитывается система из двух крыльев, в частности крыло — горизонтальное оперение. Аэродинамические силы определяются в виде матрицу коэффициентов влияния, связывающих силы в заданных узлах с амплитудами колебаний в узлах (метод конечных элементов). Каждое из крыльев рассчитывается с помощью метода потенциала скоростей как наиболее корректного из известных линейных методов. Интерференция между крыльями определяется методом потенциала ускорений и сводится к расчету возмущенных скоростей в узлах заднего крыла, индуцированных передним крылом, что для сверхзвуковых скоростей оправдано.
Замена комплексного . потенциала вещественным значительно сокращает время расчета (для 160 узлов с 25 до 5 мин на ЭВМ БЭСМ-6). Экономию времени и памяти ЭВМ дает также разбиение крыла на участки с разными типами решения: нулевое, явное, вещественное, комплексное.
Предложенная методика применена к расчету на флаттер, что иллюстрируется двумя простыми примерами, на которых показаны недопустимо большие ошибки в случае применения теории поршня при М<1,4 и хорошие результаты применения гипотезы одномерной стационарности.
1. Расчет аэродинамики крыла, колеблющегося в сверхзвуковом потоке, на практике осуществляется многими способами. Охарактеризуем некоторые из них.
Первый способ — метод потенциала скорости, когда за основу берется интеграл, связывающий потенциал скорости ср со скосом ію\
здесь приняты обозначения работы [1], которая ниже развивается на случай системы крыльев с дозвуковыми задними кромками.
Безразмерный скос w = df/dx + iShf определяется колебаниями крыла: y = bf(x, z) exp (i to t). Скос вне крыла неизвестен и определяется с помощью интегрального уравнения (1) из условия обращения
В. Г. Буньков
в нуль потенциала впереди и сбоку крыла. Сзади крыла потенциал <р ненулевой и однозначно связан со значением ср на задней кромке:
9 = 9а. к ехр [— г БЬ (д: — л3.к)]- (2)
Это соотношение следует из условия обращения в нуль давления за крылом, выражаемого с помощью формулы:
с'р — ду/дх + IБИ у, (3)
где ср—коэффициент давления, связывающий разность давлений с плотностью и скоростью: Ар = (2уигсре1Ю*.
Второй способ — потенциал на крыле по-прежнему определяется интегралом (1), но неизвестные скосы ьо в зоне влияния (впереди и сбоку крыла), а также на пелене (сзади крыла) определяются по-другому, а именно с помощью потенциала ускорений. Соответствующее интегральное уравнение, выражающее скосы через давление, имеет следующий вид [2]:
т ■
Я1
ехр (-- I Б!! хп)
Х0 ехр ( — 4>А-0/М) СОЭ/)/?
+
-l.il ' 2
М + Л
/
/ехр (— 1рЦ<и (*о ?г/Му
(И-йт.
(4)
х0 М —К
Третий способ — искомая нагрузка на крыле с' определяется через скосы ш на крыле с помощью интегрального уравнения (4). Этот способ замечателен тем, что не требует расчета зон влияния и пелены. Однако из всех способов этот способ наиболее подвержен влиянию погрешностей, поскольку интеграл (4) в отличие от (1) лишен непрерывности, а успех его использования зависит от выбора системы контрольных точек. В самом деле, если интеграл (1) дает в качестве <р непрерывную (по х, г) функцию при любом распределении скосов ш, то интеграл (4) ведет себя совершенно по-другому. Чтобы проиллюстрировать это, возьмем отдельную квадратную панель с единичным давлением на ней с'р = 1. Согласно интегралу (4) такая панель индуцирует поле скосов, стремящихся к бесконечности на боковых сторонах панели и их продолжении.
Отметим важную особенность уравнения (4). Интеграл (4) в обычном смысле при пересечении зоны интегрирования вертикалью г0 = 0 не существует. Однако в расчетах следует брать так называемое обобщенное главное значение [3]. Дело в том, что интеграл (4) получается формально как предел следующего выражения:
где у— высота, а через /(х0, z0) обозначено подынтегральное выражение (4) без знаменателя 2$. Отмеченное обстоятельство не вносит особых трудностей в процесс численного интегрирования, однако оно приводит к требованию непрерывности величины с и ее производной
на оси г0 = 0.
Предлагая первый способ в качестве одного из рабочих *, проведем некоторую его модернизацию, а именно, вводя вместо потенциала <р так называемый модифицированный потенциал Ф [4] с помощью замены
ср = Ф ехр (—¿Мрх), (6)
приходим вместо (1) к интегральному уравнению с вещественным ядром:
где W — дФ/ду(у = 0) w ехр (¿МрЧ)— модифицированный скос. Условие (2) для пелены преобразуется к следующему виду:
Ф = Ф3.кехр[г(л: — х3. к) Sh/P2]. (8)
Замена (6) приводит к вещественному виду и уравнение (4):
W
однако, к сожалению, аналог давления дФ/д^ не обладает тем замечательным свойством, как величина с' : не обращается в нуль за пределами крыла.
2, Вычисление интеграла для модифицированного потенциала скорости Ф предлагается осуществить, как и в работе ¡[1], с помощью интегрирующей сетки. Тогда интеграл (7) при достаточно мелкой сетке с нужной точностью будет равен сумме:
Ф“Е2^*Л*. (10)
/=0А=0
где
С0 = 4]/Д/3; с* = с0[(£+ 1)|/б4Г!-2£/£ + (£ — О/А^Т] , (11)
Д = 1/А^ — постоянный шаг сетки в характеристических координатах,
/у* — значение функции / =— V?соэв узлах сетки.
В отличие от [1], нумерацию узлов сетки будем связывать не с линиями Маха, а с вертикальными линиями, параллельными скорости
* Следует считать третий способ также рабочим, а сравнение их результатов — основным инструментом для их усовершенствования.
ЛдФ di
х0 cos pR R~
+ p sin pR du d% dri,
(9)
потока. Тогда из характеристической сетки ЫхЫ можно выделить 4 • N +1 вертикальных линий, содержащих различное число узлов Ь1П(т = 0,... ,4Л0, а формула (10) получает следующий вид:
™ 1ш-х
Ф=Е 02)
т~0 1=0 . ... .
где
к^2Ы~т-\-1. (13)
Формулу (12) можно рассматривать как стандартную формулу для приближенного вычисления интеграла (7). Принимаются во внимание только те узлы сетки, которые попали в область интегрирования, а для узлов, лежащих вблизи границы области; необходимо брать дробные доли от коэффициентов ¿пл. При определении дробной доли можно полагать, что каждый узел вносит вклад в интеграл не с участка площади, а с отрезка вертикальной линии, которой он принадлежит.
В расчетной модели крыла [1] предлагаются некоторые усовершенствования. По-прежнему крыло разбивается на трапециевидные и треугольные конечные элементы (КЭ). Разбиение на КЭ либо берется готовым из прочностного расчета, либо назначается. Предлагается к схеме крыла добавить боковой участок, который нужен для случая скошенной под отрицательным углом стреловидности боковой кромки (рис. 1).
Зону влияния разбиваем на три зоны: переднюю, заднюю и боковую. Эти три зоны могут присутствовать в самых разных сочетаниях, например для крыла на рис. 1 при М =1,1 имеются все три зоны. КЭ в передней и задней зонах располагаем на вертикальных полосах, общих с крылом. В боковой зоне число полос задается; а их ширину рекомендуется увеличивать по геометрической" прогрессии при удалении от крыла. По геометрической прогрессии регулируются также вертикальные размеры КЭ.
3. Для экономии времени расчета и памяти ЭВМ предварительно анализируются типы решений, ожидаемых в задаче для конкретной
формы крыла в плане и заданного числа М. Дело в том'; что в некоторых случаях решение для потенциала может быть записано в явном виде без необходимости решать систему (рис. 2,а). В других случаях решение может оказаться явным для большей части узлов и лишь на небольшой части крыла требуется решать систему уравнений сравнительно небольшого порядка (рис. 2,6, г). В свою очередь, неявная часть решения, получаемая через систему матричных уравнений, в некоторых случаях может быть чисто вещественной (рис. 2,6, в), а в других — комплексной (рис. 2,г, и) или смешанной (рис. 2Д ж, з).
Всего возможно семь типов различных комбинаций, состоящих из явного Фя, вещественного Фв и комплексного Фк решений. Необходимо также принять во внимание, что потенциал на передней и боковой кромках крыла равен нулю по условиях задачи, т. е. часть узлов имеет нулевой потенциал Фо=0.
Таким образом, решение Ф для п узлов крыла имеет следующую структуру (вектор и-го порядка):
/ Ф0\ — нулевое решение для п0 > 0 узлов,
I Ф, \ — явное решение для «я>0 узлов,
Ф ==1 Ф® | — веЩественное решение для яв> О узлов,
\ Фк /—комплексное решение для пк > 0 узлов.
При этом п = п0+пя + пв + пи — общее число узлов, образованных сетью КЭ крыла. Узлы должны быть рассортированы и пронумерованы с целью получить четыре отрезка вектора: Фо, Фя, Фв, Фк —в порядке возрастания сложности решений. Говоря о решении, мы имеем в виду матрицы, связывающие неизвестный потенциал с вектором известных скосов на крыле
Ф Я = ЛЯ1Г; Фв = Лви/; ФК = (ЛЛ. + /ЛМ)1Г, (14)
т. е. решение это не векторы Фя, Ф„, Фк, а прямоугольные матрицы ^Я> ^В> Ац, Ам.
Неизвестные скосы в зоне влияния также имеют решения разной сложности: в передних узлах передней и боковой зон скосы нулевые, а вся зона разделяется на две области с вещественным решением и комплексным (см. рис. 2):
= - Н1Г; «7, « - (Я, + Ш2) ЧР. (15)
Узлы зоны в области пелены, примыкающие к задней кромке крыла (один и тот же узел может принадлежать и крылу, и зоне влияния), имеют известные скосы, совпадающие со скосами крыла (кроме углового узла). Это условие равносильно условию нулевого давления на задней кромке, но первое удобнее для программирования, хотя при уменьшении размеров КЭ оба условия дают в пределе одинаковые решения.
Для нахождения решений Ая, Ав,..., Я2 составим систему матричных уравнений, проведя интегрирование с помощью интегрирующей сетки, накладывая ее последовательно на п^п—п0 узлов крыла и т = тя + тк узлов зоны (сетки должны быть в несколько раз мельче, чем КЭ):
Здесь указаны порядки входящих в систему прямоугольных вещественных матриц Ая, А, В, А0,..., Dy.
К этой системе добавим еще условия на задней кромке [см. уравнение (8)]:
Фс = МсФк, (17)
где Мс — комплексная матрица, состоящая из нулей и редких множителей — не более одного в строке (индекс «с» означает след) следующего вида:
pir = cos я —f— ¿ sin a; a = sh (Xir — A'3 Д/р2,
где і — номер цепочки узлов в пелене, выходящей из і-го узла на задней кромке, а г — номер узла в цепочке; Xir—X3i — расстояние от задней кромки до рассматриваемого узла. Для удобства программирования рекомендуется нумеровать узлы крыла так, чтобы опорные узлы на задней кромке, из которых выходят цепочки узлов пелены (в сумме насчитывающих т3 штук) с ненулевыми потенциалами, оказались в конце вектора W, а упомянутые т3 узлов пелены — в начале отрезка вектора WK из тк узлов. Все умножения на матрицу Мс, которые встретятся ниже, следует заменить простым умножением строк на множители Ціг.
Система (16) — (17) решается так:
ФК = Л0^+ Яо^в + Я.^к, Со Ш + Я0 + Я, 1ГК = Жс Фк,
откуда с учетом решения (18) следует:
где
где
Ф , = Аг1Г+В>\Ук, СгХУ + О^-М, Ф„ Л,=Л0-Я0Я; -Сх = С0 — О0Н\
С2 \Г + 021ГК = 0,
С2 — Сі — Л^Л); = — МсіЗі!
Шк = ~Ої1С2№= — (Н1 + ¿//*) ИГ.
В результате получаем значения #1 и Я2.
В операции (19) в отличие от (18) появляются действия с ком-плекными матрицами, однако и здесь эти действия следует привести к аналогичны^ для вещественных матриц; в частности, подметив особенность матрицы 1>2, заключающуюся в том, что только первые т3 ее строк содержат мнимые части, можно ускорить два из трех перемножений в отношении т3\тк (обращение комплексной матрицы сводится к трем перемножениям и двум обращениям).
И, наконец, получаем решение для величины Фк [см. (14)]:
Получив решение для модифицированного потенциала, т. е. комплексную матрицу влияния Л(ф = Л№), легко переходим к обычному потенциалу: <р — Ра>, для чего каждый из коэффициентов а,к матрицы А умножаем на фазовый множитель согласно замене (6) —(7):
Так получается матрица Р, выражающая потенциалы ф в п узлах, образованных сетью КЭ, через возмущения скорости в них и).
При переходе от потенциалов ср в узлах к нагрузкам Р в узлах (Р = = <?ау) применяем принцип, описанный ранее [1], но обращаем внимание на то, что составлять матрицу перехода (2 от потенциала ф к силам Р(Р=фф) не следует, чтобы потом не тратить время на перемножение (й = С}Р). Вместо этого следует применять операции типа линейных преобразований строк: каждая строка й получается линейной комбинацией нескольких строк Р в соответствии с нумерацией узлов, окружающих рассматриваемый. Чаще всего в образовании строки О участвуют девять строк Р (к узлу обычно примыкают четыре КЭ).
4. Расчет целого летательного аппарата построим на комбинации двух несущих поверхностей (рис. 3). Как правило, передняя поверхность — это крыло, а задняя — стабилизатор. Однако возможно и об-
/>к = ехр [¿Мр (хк- х,)\, (///у,;/ = Р),
(21)
2—«Ученые записки» № 3
17
ратное (рис. 3,д), поэтому будем называть переднюю поверхность «крыло 1», а заднюю — «крыло 2». Основной рабочей гипотезой является предположение, что крыло 2 не влияет на крыло 1. Действительно, для очень многих схем ЛА даже при М=1,05 (предельно низкое число М, при котором еще можно надеяться на пригодность применяемой линейной теории) линии Маха достаточно круты (1§ф = 0,32), чтобы, будучи
проведенными от передних углов крыла 2, не задеть крыла 1. Исключение составляет схема ЛА с поворотным крылом, когда крыло подходит близко к стабилизатору л (рис. 3,в). Граница между крыльями / и 2 должна выбираться из условия отсутствия влияния второго крыла на первое, как это показано на схемах гиена рис. 3.
Имея решение для каждого крыла в отдельности:
Рг = 02Щ, (22)
и рассчитав скосы шИНд, индуцированные крылом 1 в узлах крыла 2, можем написать решение для системы:
Р2 = е2(и>2 — о>инд). (23)
Индуцированные скосы будем рассчитывать с помощью формулы (4), которую представим в следующем виде:
ву = -^- ¿ = е-1&Ъх<1{е ™ соъ рР—]. (24)
^ I/ */ К V Х0 )
\ '
Замена формулы (24) на (9) никакого выигрыша не даст, так как здесь речь идет о простой квадратуре.
Учитывая плавный характер функции g (эта функция непрерывна вместе с ее производной) и ее монотонность для достаточно малых областей интегрирования, какими являются К.Э, можем написать:
„lÜÜlifÄiü,, (25)
ъ J J Zq R
где x*, z* — некоторая точка внутри КЭ, наиболее правдоподобное положение которой находится в геометрическом центре тяжести КЭ (к такому выводу приводит замена функции g(g, rj) ее линейным приближением, а остальной части подынтегрального выражения (24) — усредненной по площади КЭ постоянной).
В интеграле (25) мы допускаем еще одно приближение, заменив функцию с' константой, равной суммарной силе, действующей на КЭ, деленной на площадь КЭ. Конечно, точность вычисления интеграла (24) была бы выше, если бы этот интеграл удалось взять в явном виде для линейной интерполяции функций ср и g, или их произведения, но, к сожалению, этого сделать не удается.
Таким образом, для интеграла на КЭ (или на фрагменте КЭ, отсекаемом конусом Маха) имеет формулу:
w = g'(xa.r,znT)cpI; g' = — JL; /= — ff X-d2^r‘ . (26)
sK3 го x
Здесь комплексный множитель g характеризует запаздывание по фазе и при нулевом числе Струхаля превращается в единицу. Интеграл I, характеризующий стационарную задачу, берется в явном виде:
/ = ^ + ,п7Т^Т/ + т?; т=/1Р2-^1. (27)
где cp- arcsin ^-^1... или In+ % + (PS И),
t — тангенс угла стреловидности передней (задней) кромки КЭ, h — отрезок на оси х0 при пересечении упомянутой кромкой.
Матрицу усредненных давлений с’р в КЭ выражаем через матрицу потенциалов. В конечном итоге получаем матрицу индуцированных скосов в заданной сети контрольных точек (назовем ее вспомогательной
сетью узлов — ВСУ), как функцию скосов ш крыла. Узлы ВСУ нельзя
задавать произвольно, их следует задавать на средних линиях трапеций — КЭ. Пересчет от узлов ВСУ к узлам крыла 2 делаем с помощью линейной интерполяции. В результате вектор шИН1Т определится некоторой прямоугольной матрицей В\
™иНд =Bwl. (28)
Описанная схема расчета ЛА может быть обобщена на большее, чем две, число несущих поверхностей и к тому же расположенных не в одной плоскости, а произвольно в пространстве, для чего формулу (4) следует заменить другой, определяющей пространственный скос потока от заданных давлений. Расчет же каждой отдельной поверхности следует проводить методом потенциала скорости.
5. Методика применена к расчету симметричных форм флаттера самолета (или антисимметричных, но без учета аэродинамики киля). При описании органов управления (элеронов, рулей) сеть КЭ строится
таким образом, чтобы границы рулей проходили по возможности посредине КЭ — это связано с тем, что скосы потока w задаются в узлах КЭ, но скосы терпят разрыв на границе рулей. Упомянутый прием устраняет неопределенность и заменяет разрыв плавным переходом.
Проиллюстрируем методику на двух простых примерах: флаттер квадратной и трапециевидной пластинок. Попутно сравним полученные результаты с теорией поршня и гипотезой одномерной стационарности (ГОС), которые часто применяются на практике. По теории поршня (вернее, по ее линейному приближению) давление в любой точке крыла равно:
р(х, z) = — с’р(х, г)2?ъ*(ду1дх-\-у№)\ ср (х, z) = 1/М, (29)
где у — деформация в рассматриваемой точке.
Основное свойство этой формулы — ее одномерность — используется в способе, который назовем гипотезой одномерной стационарности. Будем определять давление по-прежнему формулой (29), но' коэффициент с' (х, z) зададим как функцию точки х, г, исходя из стационарного решения или экспериментального распределения давлений. Поскольку речь идет о задачах аэроупругости, т. е. в основном об устойчивости, то правильнее в качестве величины ср(х, г) брать не само давление, а приращение давления при приращении угла атаки, т. е. вариацию давления. Такое определение с'р приводит к возможности использования ГОС для толстых тел под большими углами атаки (нелинейность) с учетом скачков и т. п.
Допустимость применения ГОС объясняется тем, что роль нестационарное™ для сверхзвуковых крыльев не особенно важна; так, учет чисел Струхаля (реальные числа Sh при флаттере обычно равны 1—2) дает поправку обычно около 10% . В то же время ГОС позволяет более правильно задать распределение давления по сравнению с расчетом по линейной теории, в частности это касается органов управления. То, что в ГОС берется давление для жесткого крыла, а решается задача об упругом крыле, не имеет решающего значения, так как деформации обычно плавные.
Первый пример расчета. Квадратное консольное крыло — пластина постоянной толщины. Задаем на нем сеть из т= 10X10 КЭ или 11X11 = 121 узла (рис. 4), M = y^2, Sh = 0. Точное решение для ср равно единице в треугольнике ABC и агссоэ ф/я в треугольнике BCD,
где "ф линейно меняется слева направо от ф = —1 на линии ВС до ф=1 на линии BD. Расчет давлений в узлах дает ошибку в среднем около 3% (величины а, су, Хф определяются с погрешностью 1%).
Теперь представим это крыло как систему двух крыльев: переднее 1 и заднее 2 (см. рис. 4) Расчет для крыла 2 (на кры-
ле 1 ничего не меняется) дает на этот раз
2 по-прежнему в большинстве узлов ошибку 3%, однако 10 узлов имеют ошибку около 10,% (ошибка с”, Хф системы составляет 2%).
Рассчитаем квадратное крыло на флаттер, для чего назначим конкретные Рис. 4 размеры, массу и жесткость. Сторона квад-
А
* щ
/ '1
/ Л 7 1
/ / / / I
/ / t / /
/ / / / / /
/ t / ' /
/ J / /
с / t t
рата равна 10 м. Крыло — однородная трехслойная панель (обшивка и заполнитель). Толщина обшивки 6 = 0,2 см, ее модуль и плотность: Е = 2- 107 Н/см2; р = 25 г/см3 (в этой плотности учтена еще и масса заполнителя). Плотность воздуха р0 = 1,25 • 10-3 кг/м3, иэв = 300 м/с.
Приведем результаты расчета. Частоты в пустоте равны оэ1 = 12,9, (02 = 31,1 с-1, частота флаттера ЮфЛя=20 с-1 (БЬ^ОД). Критическая относительная плотность, сосчитанная по разным теориям, показана в таблице:
Одно крыло Система крыльев ГОС Теория поршня Расчет [5]
Sh = 0,1 1,53 1.62 (6%) 1,33 (-10%) 8,3 1,54
Sh = 0,44 1,47 1.41 (-4%) (число КЭ равно 36)
Sh = 3 1,36 ! ,44 (6%)
Время ЭВМ 4'40" 5'30" 18" 18" 40"
Сравнение результатов говорит о правильности предлагаемой методики расчета системы крыльев (ошибка 4—6%). ГОС дает ошибку 10% (треугольник BCD на рис. 4 разбивался на пять участков с постоянным значением с'—это участки с границами г|5= ± 1; ±0,6, ±0,2).
Теория поршня приводит к ошибке в пять раз из-за неправильного положения фокуса. Выигрыш времени дает расчет по методике [5], однако она принципиально непригодна для дозвуковых задних кромок, так как интеграл (4) вычисляется там по формуле Гаусса и имеет конечную величину только спереди и сбоку от крыла.
Расчет трапециевидного крыла на флаттер. Рассмотрим крыло-пластину из работы |[6]. Кратко запишем параметры крыла: Ь = 200 см, /=120см, tgx=l,25, tg%i = 0, h = 2 см, 6 = 0,05 см, Е = 7 ■ 10е Я/см2, р = = 4,8 г/см3, частоты coi = 122; 361 с-1, число М—1,41. Разбивка крыла на конечные элементы сделана двумя способами: косоугольным
(рис. 5, а) и прямоугольным (рис. 5.6). Получившееся при этом число
КЭ такое: 120 и 128 (число узлов соответственно 143 и 156 — порядок матриц определяется числом узлов). На рис. 5 показано, как представляется данное крыло в виде системы двух крыльев.
В таблице графы К1 и К2 означают косоугольную разбивку в виде одного крыла и в виде системы двух крыльев. Аналогично П1 и П2 — для прямоугольной разбивки.
К1 К2 П1 П2 ГОС Поршень 15]
Sh =0 0,78 0,62 __ 0,44 3,32 0,80
Sh = 1 0,79 0,62 0,79 0,76 — — (число КЭ равно 35)
шфл 250 220 250 260 270 440
МБ 82 50 123 58 13 13
Время 3'26" 3'05" 4'48" 3'35" О'ЗО" О'ЗО" 0'4С"
Как видно, расчет в виде одного крыла дает устойчивый результат для любой разбивки и любого числа Струхаля. Ошибка для системы двух крыльев меняется от —18 до —4% в зависимости от способа разбивки. Теория поршня дает по-прежнему огромную ошибку — в четыре раза. Однако и ГОС на этот раз дает ощутимую ошибку 44% (в пересчете от критической скорости флаттера это означает ошибку 20%).
ЛИТЕРАТУРА
1. Буньков В. Г. Расчет методом конечных элементов аэродинамических сил на колеблющемся крыле в сверхзвуковом потоке. — Ученые записки ЦАГИ, 1981, т. XII, № 4.
2. Watkins С. Е., Berma д I. Н. On the kernel function of the integral equation relating lift and downwash distributions of oscillating wings in supersonic flow. NACA Report 1257, 1956.
3. Общая теория аэродинамики больших скоростей/Под ред. проф.
Сирса. — М.: Воениздат, 1962.
4. Джонс (W. P. Jones), Anna (Kari Арра). Применение метода градиента потенциала к нестационарной сверхзвуковой аэродинамике. —
РТК, AIAA, vol. 15, 1977.
5. Кузьмина С. И. Расчет нестационарных аэродинамических характеристик тонких крыльев конечного удлинения при сверхзвуковых скоростях полета. —Труды ЦАГИ, 1981, вып. 2118.
6. Буньков В. Г. Расчет оптимальных флаттерных характеристик методом градиента. — Труды ЦАГИ, 1969, вып. 1166.
Рукопись поступила 1/VU 1982 г.
