Системный подход к расчету оптимальной формы низколетящего крыла методом релея-ритца

Задачи поиска оптимальных решений при проектировании сложных технических систем на основе математического и численного моделирования или различных эвристических подходов возникают в различных областях человеческой деятельности. Сюда можно отнести и проблему оптимального проектирования летательных аппаратов (ЛА) различных компоновочных схем. Актуальность подобных задач была озвучена еще в начале XX века на заре авиастроения. Не менее важным направлением поиска оптимальных аэродинамических форм является проблема устойчивого движения ЛА. Особенно для низколетящих ЛА, движущихся вблизи границы раздела сред (ЛА на динамической воздушной подушке). Всесторонними исследованиями этих задач в 70-80-х гг. XX в. занимались представители школы А. Н. Панченкова [1]. В отличие от наиболее распространенного подхода, в котором для обеспечения устойчивости низколетящего ЛА применялась самолетная схема со стабилизатором, размещенным вне зоны влияния твердой границы (см. напр. [2]), А. Н. Панчен-ков и его коллеги занимались разработкой аппаратов по схеме «утка». Устойчивости ЛА самолетной схемы добиваются выносом хвостового ста -билизатора из зоны влияния опорной поверхности вверх [2]. У ЛА, построенных по схеме «утка», той же цели стремятся достичь совместной работой носового стабилизатора и основного крыла у опорной поверхности. Крыло и стабилизатор, связанные между собой, действуют в этом случае как единая система, подверженная влиянию границы раздела сред. Если в первом варианте компоновочной схемы влиянию твердой границы подвергается единственное крыло, то во втором - пара связанных между собой несущих поверхностей (крыльев). В работе [3] показано, что для подобных систем при исследовании критерия самостабилизированности целесообразно проводить аналогию с остойчивостью корабля. При этом главной характеристикой стабилизированности будет метацентрическая высота (так называемый «корабельный» подход).

Решить проблему обеспечения стабилизации летательных аппаратов вблизи опорной поверхности возможно двумя способами:

• -    разработать специальную автоматическую систему стабилизации;

• -    разработать специальную геометрию форму несущих поверхностей летательного аппарата, а также их компоновки в аэродинамической схеме. Такой способ называется аэродинамической стабилизацией.

Для малых и средних летательных аппаратов проблема самостабили-зации может быть решена выбором специальной конфигурации крыла.

На больших летательных аппаратах выигрыш от применения аэродинамической стабилизации заметно снижается, оснащать их системой автоматической стабилизации предпочтительнее. Также хорошие перспективы у комбинированных систем [1].

В случае использования аэродинамической стабилизации либо при разработке комбинированных систем стабилизации необходимо выбрать наиболее подходящую (оптимальную) конфигурации летательного аппарата и его главных несущих элементов - крыльев.

Исследуемая система состоит из низколетящего крыла конечного размаха (несущая поверхность), обтекающего потока идеальной жидкости (без учета вязкости) и твердой поверхности (это может быть грунт, лед и т. д.) (рис. 1). Полагаем, что несущая поверхность находится на сверхмалом отстоянии от границы раздела сред различной плотности, составляющем менее 0,1 хорды крыла.

В рассматриваемой оптимизируемой системе (крыло, твердая граница, поток) объектом проектирования является несущая поверхность, так как на стационарный поток и границу раздела сред повлиять не представляется возможным.

Рис. 1. Изолированное низколетящее крыло над границей раздела сред.

Известно, что аэродинамика крыла на сверхмалых отстояниях от твердой границы определяется формой его нижней поверхности. В связи с этим главную роль в формировании АДХ низколетящего крыла играют геометрические характеристики именно нижней поверхности [4].

Эффект влияния твердой границы определяется относительной высотой полета (отстоянием от границы), отнесенной к хорде крыла. При движении вблизи твердой границы существенно растет подъемная сила крыла и снижается его индуктивное сопротивление. Если высота полета меньше 0,1 длины хорды несущей поверхности (сверхмалое отстояние), влияние границы раздела сред настолько возрастает, что аэродинамическое качество увеличивается до 40-50.

Крыло непосредственно взаимодействует со стационарным потоком, в результате чего в потоке формируются возмущения. Параметром возму щения служит отношение е =

a

, где a

— местный угол атаки,

h — местное относительное отстояние в точке g , лежащей на поверхности крыла.

В зависимости от величины параметра е в динамике низколетящего крыла выделяют следующие варианты теорий:

• 1)    е <<  1 — теория малых возмущений;

• 2)    е ~ 1 — теория умеренных возмущений;

• 3)    е >>  1 — теория сильных возмущений.

В рассматриваемой системе возмущения полагаются умеренными.

Одним из первых результатов в области оптимизации крыла стала теорема Мунка о крыле в безграничной жидкости с минимальным индуктивным сопротивлением при заданной подъемной силе. Было показано, что такое крыло имеет эллиптическую по размаху циркуляцию.

Первые результаты по определению геометрии плоского крыла вблизи границы раздела сред с минимальным индуктивным сопротивлением получил П. Халлер в 1936 году: для плоского низколетящего крыла минимум индуктивного сопротивления достигается при параболическом по размаху распределении циркуляции.

Для низколетящих крыльев в 80-х гг. прошлого века была разработана теория оптимальной несущей поверхности (ТОНП) [4]. Методы ТОНП опираются на теорию потенциала ускорений, квадрупольную теорию крыла, асимптотическое программирование, результаты исследований экстремальных задач теории крыла в ограниченной жидкости.

ТОНП характеризуется тем, что большинство ее задач относится к классу некорректных экстремальных задач вариационного исчисления. Некорректность вызвана тем, что количество граничных условий, которые накладываются на экстремаль, превышает порядок дифференциального уравнения необходимого условия экстремума. Это их специфическая особенность.

Второй особенностью служит их сложность, что затрудняет получение аналитических решений. В [5] был предложен подход к получению таких решений, но в общем случае целесообразно применять численные методы оптимизации. В работе для этого используется метод Ритца, который представляет собой частный случай метода пробных функций. Согласно ему, задача по определению минимума соответствующей функции преобразуется в задачу нахождения корней системы алгебраических линейных уравнений. Численное решение такой системы не вызывает затруднений даже при числе параметров n ~ 100 - 200.

Постановка задачи. Одной из принципиально важных проблем проектирования ЛА на динамической воздушной подушке является обеспечение его самостабилизации. Для малых и средних ЛА решением проблемы самостабилизации может быть выбор специальной конфигурации крыла. При этом одним из критериев самостабилизации судна на динамической воздушной подушке является метацентрическая высота [3]:

H M

C9 Г с9 y Cm с h ) _m с h

Cy J

С  С9

C y V C y

Чем она больше, тем лучше.

Здесь C_y — коэффициент подъемной силы крыла, C 9 — производная

Cy по углу тангажа, Cy — производная Cy по отстоянию, Cm — коэф- фициент момента тангажа, С^ — производная Cm по углу тангажа, Cmh — производная Cm по отстоянию.

К сожалению, в настоящее время не представляется возможным решить эту задачу комплексно с учетом всех аэродинамических производных, поэтому ограничимся частным случаем задачи: h у h ^ min, Cy поскольку по физическим соображениям можно предположить, что вблизи границы раздела сред именно производные по отстоянию будут играть главенствующую роль. Она формулируется следующим образом: определить геометрическую форму крыла с максимальным значением коэффициента аэродинамической производной подъемной силы Cy по отстоянию h , если значение коэффициента аэродинамической производной мо-

C m = const;

C y ^ max.

мента тангажа C_m по отстоянию h фиксировано:

При решении такого рода изопериметрических задач в качестве функционала используется сумма: Ф ₁ + рФ 2 . В нашем случае Ф ₁ = C y , Ф ₂ = C^h_m , р - произвольная константа.

Задача Cy ^ max в приближении умеренных возмущений была решена в [5]. В этой же работе для умеренных возмущений получен коэффи-h циент производной момента тангажа по отстоянию Cm :

	Ch = лщ г F ( y ) dy _ лтц r F i ( y ) d y + m 2 H o H 0 0 H y ) y’ 8 Л 2 0 Н ,2 J H ( y ) +s OH i « y J h ( y F ( y i ) F y i ( y ) dy i - 8 H H 0 0 л ( y ) y 8 л2 H 2 ii a ‘Х ( * ) F 2 ( ’ ) 8 Z o H о о о a 2 ( X ) f F 2 ( y )       i )!,, i                              Xdxdy , H ( H 2 ( y ) H ( y H i H i2 J
При этом:	C y =      i F ( ’ ) dy        i 2 [ a iX ( X ) H ( ’ ) F 2 ( ’ )+ 2 H 00          8 H 000 . i —2( \ Fy ( y )      i       J Я + H "  'X\ H ( y )   \ H 2 J] dxdy’

1 y 1 dy Т7М-Г У 1 dy 1 ^де F ⁽ y’^H ( У 1 )’ F ^{1 ( y} ) = H ( У 1У

Здесь X — удлинение крыла, X ( y ) = X ( y ) / X — безразмерное местное удлинение крыла ( X ( y ) = b / a ( y ), b — полуразмах, a ( y ) — местная полухорда крыла в плане), X ₀ = b / a ( 0) — относительное удлинение крыла, a i — полный угол индуктивного скоса потока по задней кромке, a i ( x ) = a i ( x ) / a i — местный безразмерный угол индуктивного скоса потока, H - относительное по размаху безразмерное отстояние несущей поверхности от границы раздела сред, H ( y ) = H (0, y )/ H о — форма задней кромки крыла, H 1 = H (1) a x ( x ) — производная от a i ( x ) по x , F_y ( y ) — производная от F ( y ) по y , F₁ y ( y ) — производная от F₁ ( y ) по 2

y , X i = X (1), n = J a / ( x i ) dx . о

Следует отметить, что данные функционалы получены для т.н. «крыльев Мунка», характеризующихся тем, что их местные углы индуктивного скоса потока постоянны по размаху при движении вдоль вихревых линий, подобных по форме в плане входной кромке крыла (заднюю кромку полагаем прямолинейной в плане). Этим обеспечивается поиск решения на классе крыльев с минимальным индуктивным сопротивлением при заданной подъемной силе [5]. Таким образом, задачу определения геометрии крыла с максимальным значением коэффициента аэродинамической производной подъемной силы Cy по отстоянию h при фиксированном значении коэффициента аэродинамической производной момента тангажа Cm по отстоянию h решаем именно на крыльях Мунка. В общем случае оптимизация выполняется по параметрам a/x, X (y) и H(y) тематическая постановка задачи имеет следующий вид:

, а ма-

^ = max {Ф1 + рФ 2};

W

W1 = {a/ (x),X (y), H (y)};

a(0) = 1, ax (0) = 0, a(2) = 0;

X (1) = kX(0);H(0) = 1,H(1) = %0;

J =1. H (y ) dy=e

Здесь Ф ₁ = C h , Ф ₂ = C^ m , p — произвольная константа.

Благодаря специальному выбору параметра a i ( x ) экстремальная задача (6) станет эквивалентна проблеме экстремизации вклада нелинейного слагаемого в исходный функционал (7). Как указывалось ранее, задача является некорректной. Третьим («лишним») граничным условием в ней является постулат Жуковского-Чаплыгина u ( 0 ) = 0.

Граничные значения для параметра a -( x ) , которые входят в эту постановку, вытекают из условий, наложенных на погонную циркуляцию Г ₀ ( x , у ) : не равная нулю нагрузка на крыле, постулат Жуковского-Чаплыгина, а также условие безударного входа – соответственно. Изопериметрические условия для параметров X ( у ) и H ( у ) следуют из ограничения на величину подкупольного объема V ( Q ) = const .

Как было указано, геометрические характеристики крыла складываются из характеристик крыла в плане, характеристик профиля и характеристик задней кромки. Поэтому задачу определения оптимальной геометрии крыла целесообразно разбить на три частные подзадачи:

• 1)    задача определения оптимального распределения индуктивных уг-

• лов скоса потока по хорде;

• 2)    задача определения оптимальной геометрической формы крыла в

плане;

• 3)    задача определения оптимальной геометрической формы задней кромки крыла (вид спереди).

Рассмотрим задачу определения оптимального распределения углов скоса потока по оси 0 x (по хорде) для несущей поверхности крыла с максимальным значением коэффициента аэродинамической производной подъемной силы C y по отстоянию h , когда значение коэффициента аэродинамической производной момента тангажа C m по отстоянию h фиксировано. В формализованном виде постановка исследуемой оптимизационной задачи такова:

^1 = minj(u2 (t)(1 + 2pt) + u2 (t)(P - 2ptp))dt;

u ( 0 ) = 1, u ( 0 ) = 1, u ( 1 ) = 0.

Здесь u = a i ( t ), u ( t ) — производная u ( t ) по t ,

Р =

4 1 ( F y C y ) ;^{2 j л} о

l ^{л (} У ) Л Hl

dy

j F ² ( у ) ! ( у ) dy

;

r x ( F y ⁽ y ) ^F1 y ( y ) ² = о ^л p =—

l  л2 (у)    л(у Л H j F (у)F1 y (у) dy 0

—

^

;

1 7

dy

;

p — произвольная константа,

Вместо переменной x используется формальная переменная t = ^ХД . Производная по переменной t обозначается точкой.

Применим численный способ нахождения минимума функционала.

Общая схема численного решения оптимизационных задач заключается в сведении задачи

^{ф [} y ^{( x )]} = ^{inf q [} у ^{( Х )]} ; у ^{( Х} ) ^{у ( Х ) е} Y к поиску минимума функции многих переменных [5].

В исследуемом случае краевые условия задачи оптимизации неоднородны.

Для решения оптимизационной задачи (6) будет применен численный метод Ритца, также называемый методом пробных функций.

Описание метода Ритца. Используемый здесь метод основан на редукции исходной экстремальной задачи к последовательности конечномерных задач, решения которых сходятся к искомому. С этой целью, в соответствии с методом, строится соответствующая минимизирующая последовательность.

В качестве примера можно рассмотреть задачу определения минимального значения функционала V [у] в классе функций M. Сделаем предположение, что конечный инфинум ц значений функционала существует, а также, что в классе допустимых функций должны существовать такие функции, на которых значения исходного функционала конечны. Отсюда можно сделать вывод, что должна существовать такая минимизирующая последовательность функций у1,у2,...,yn,..., которая удовлетворяет условию lim V[ Уп ] = ц.

n ^^

Если последовательность у 1 ,у ₂ ,..., у_п ,... имеет предел у *, то это и определит решение задачи, ибо выполняется предел:

V[У ]= lim V[Уп ]• п ^ю

Для получения такой последовательности следует выбрать базис в виде системы функций ф ₁, ф 2 ,..., ф _п ,..., которые должны удовлетворять условиям :

• -    принадлежать классу M ;

• -    их всякая конечная линейная комбинация у_п = с₁ ф ₁ + c₂ ф 2 + ... + с_п ф _п также принадлежит классу M .

Чтобы минимизировать функционал

V [ z ] = J[ Р (t) y 2 (t) + k (t) y (t) +

на множестве функций M = { y ( t ) y ( t ) g C [0 ,1 ] , У ( 0 ) = y ( 1 ) = 0 } следует сфор- n мировать конечную систему таких линейно-независимых функций { ф , } i = ₁, чтобы они удовлетворяли краевым условиям ф _/ ( о ) = ф _/ ( 1 ) = 0, i = 1, n . Это позволяет заменить исходную задачу на более простую, связанную с минимизацией функционала на совокупности линейных комбинаций:

n

ф(t )=Е ciФi(t )

i = 1

Решение задачи определения оптимальной формы несущей поверхности методом Ритца. Перейдем к задаче минимизации рассматриваемого функционала на линейных комбинациях, для чего подставим ф ( t ) в функционал:

n

n

n

n

n д V

Необходимое условие экстремума имеет вид--- = 0.

5 C j

Продифференцировав функционал V по параметру c _j , получаем:

дV 1 I

"Т" = Н2Р(t)Ёc^^- (t)фj (t)+ k(t)фJ (t)+ д ci al           i=1

' j 0 I i = 1

n

+ 2q(t)^c^i (t)ф} (t) + f (t)ф} (t)^dt = 0.

i = 1

После преобразования необходимое условие экстремума примет вид:

1/

Ё j (2 p (t ^i (t ^j (t)+2 q (t р (t )v j (t ))dt ci + i=1

1                                    1

+ j f ( t V j ( t d + j k ( t V j ( t "W - 0                   0

Запишем полученную систему в матричной форме: A • c = b .

Здесь

^/„/х aj =№ p(t pi(t pj(t)+2 q(t Mt pj(t M bj = -j(f (t >j (t)+k (t pj (t ydt .

0                                                          0

Ниже описаны этапы, которые необходимо выполнить для нахождения оптимальной геометрической формы несущего крыла методом Ритца:

• -    приведение неоднородной оптимизационной задачи (6) к эквивалентной ей задаче с однородными краевыми условиями;

• -    выбор базисных функций, в качестве которых можно использовать кусочно-линейный или кубический интерполяционный базис;

• -    построение экстремизирующей последовательности;

• -    восстановление геометрической формы несущей поверхности.

Результаты применения метода Ритца к поставленной задаче оптимизации. В качестве базисных функций будет использован кусочнолинейный базис. Далее необходимо построить систему кусочно-линейных функций, для чего отрезок [0,1] будет разбит на n отрезков

0 = 10 < 11 < ••• < tn < tn+1 = 1 длиной hi = t,+1 -1,, i = 1,n . Базисные функ- ции записываются следующим образом:

0,             если 0 <  t <  ti.

P i ( t ) =■

1 /       x , если t i - 1 <  t <  t i hi- \ ( t - t i - 1 ) J- ( t l + l - t ) .    если t > < t <  t i + V h - 0,              если t i + 1 <  t <  1

Будем считать, что длина отрезков h i = h = 0,1, t i = 0,1 i , i = 0,9 . Известно, что i -я функция отлична от нуля только на промежутке [ t i - 1 , t i _{+ 1} ] .

Поэтому справедливы равенства:

P i ⁽ t ) Р j ⁽ t ) = ° и Р i ⁽ t ) Р j ⁽ t ) = 0, i , j = 1, n , j Ф i - 1, j Ф i , j Ф i + 1.

Следовательно, матрица системы уравнений, содержащая неизвестную величину c_i , при таком выборе базисных функций будет трехдиагональной. Ненулевые значения будут иметь элементы:

aii = I4,i + I4,i+1 + I2ii + I3,i, i = 2, n, ai, i+1 = - I 4, i+1 + 11, i, i = 1, n - 1, ai,i-1 = -I4,i + 11,i-1, i = 2, n, bi = 15,i + 16,i + 17,i + 18,i, i = 1, n, где I - интегралы вида

1 м =

^- t! \t ^- t i ) q ⁽ t ) dt .

Эти интегралы необходимо вычислить для построения трехдиагональной матрицы системы уравнений. Далее, на основе полученной системы линейных уравнений вычисляются значения пробных функций и форми руется искомая оптимальная геометрия крыла.

Описанный выше численный алгоритм был использован для решения оптимизационной задачи определения геометрии крыла с максимальной подъемной силой [7] с целью его верификации. Аналитическое решение задачи определения геометрии крыла с максимальной подъемной силой было получено ранее в работе [4]. Численный алгоритм был реализован в компьютерной программе. Качественное совпадение результатов подтвердило работоспособность алгоритма и достоверность полученных с его помощью данных.

Найденная форма экстремальной функции позволяет определить искомую оптимальную геометрическую форму поверхности крыла. Угол атаки определяется по формуле [4]:

a x , У ) = ^a i ( x )^{- Л} 2 ( У ) H ( У ) ^F ( У ) d a ^{i (} x ) dx

.

Геометрическая форма крыла определяется по формуле [4]:

П(x, У) = Ja(x)dx - Л2(у)Н(у)F(у) da(x). о                            dx

Приведем исходные данные для выполнения расчетов с помощью компьютерной программы:

р=Р=п%, i~=1, л(у)=з.7, f(у)=^-т; р£[2;8].

Геометрическая форма несущей поверхности, которая была получена с использованием численных методов, требует наложения на радиус кривизны дополнительного изопериметрического условия.

Формы нижней поверхности крыла, максимизирующего коэффициент аэродинамической производной подъемной силы по отстоянию при фиксированном значении коэффициента аэродинамической производной момента тангажа по отстоянию при коэффициенте р = 2 и р = 8 в сечениях у = 0, у = 0.5 и у = 1, полученные методом пробных функций, показаны на рис. 2.

Полученная форма поверхности качественно повторяет форму поверх -ности крыла с максимальным градиентом подъемной силы по отстоянию [4].

Рис. 2. Форма нижней поверхности крыла в сечении у = 0 (график 1), у = 0,5 (график 2), у = 1 (график 3).

Заключение

О работоспособности предлагаемого алгоритма свидетельствует тот факт, что для верификации он был использован при определении формы несущей поверхности крыла Мунка с максимальным коэффициентом подъемной силы. Эта же задача была решена ранее в работе [4].

В данной задаче экстремаль и геометрическая форма несущей поверхности летательного аппарата, полученные с использованием численного алгоритма на основе метода Ритца, оказались качественно подобны экстремали и геометрической форме нижней поверхности крыла с максимальным градиентом подъемной силы по отстоянию, полученной в этой работе. Следовательно, предлагаемый алгоритм может быть использован для оптимального проектирования внешней геометрии несущих комплексов летательных аппаратов, движущихся вблизи границы раздела сред, и позволяет находить решения для задач оптимизации, аналитическое решение которых невозможно либо связано с большими трудностями .

Решение рассматриваемой в статье изопериметрической задачи позволило восстановить форму крыла, идентичную форме, полученной из решения задачи максимизации градиента подъемной силы по отстоянию h [4].