METHODOLOGY FOR ASSESSMENT OF THE WEIGHT CHARACTERISTICS OF THE CARGO AIRCRAFT WINGS USING MULTIDISCIPLINARY NUMERICAL SIMULATION AND ANALYTICAL METHODS OF MULTI OBJECTIVE OPTIMIZATION

Full Text

ВВЕДЕНИЕ Задача выбора основных геометрических параметров схемы самолёта на начальных стадиях проектирования формулируется в терминах нелинейного математического программирования [1, 2, 3, 4], что предполагает минимизацию (максимизацию) комплексной целевой функции-свёртки при удовлетворении геометрическим и функциональным ограничениям. Часто в качестве таких свёрток принимают взлётную массу самолёта, коэффициент топливной эффективности [1, 3] или более общие критерии, как себестоимость перевозок или прибыль от транспортной операции [5]. Важнейшим слагаемым таких целевых функций является масса конструкции планера самолёта, оказывающая непосредственное влияние на её итоговое значение. В связи с этим точность оценки массы конструкции планера самолёта, а в частности, крыла, несёт ключевое значение для достоверного определения основных лётно-технических и технико-экономических показателей будущего летательного аппарата (ЛА) ещё на начальных стадиях проектирования. Традиционные методики проектирования самолётов, сложившиеся в ведущих конструкторских бюро за длительное время, используют так называемые весовые формулы для оценки массы конструкции планера самолёта, которые базируются на результатах обработки и обобщения обширного количества статистического материала в совокупности с простыми аналитическими соотношениями сопротивления материалов, связывающими вес материала элементов конструкции и величину внутреннего силового фактора при типовых видах силовой работы. Качественная систематизация статистических данных и их скрупулёзный анализ позволили выстроить надёжный инструмент расчёта массы существующих конструкций с возможностью отслеживания влияния на её величину различных параметров ЛА, его особенностей, применяемых материалов, условий эксплуатации, предъявляемых требований и прочих факторов. Работа с такими математическими моделями требует минимальных временных и вычислительных затрат. По сути это чётко определённые аналитические зависимости, к которым могут быть применены аналитические непоисковые методы нелинейного программирования для использования их в задачах оптимизации параметров облика самолёта. Разработкой весовых формул занимались практически во всех ведущих инженерных коллективах по всему миру. Подробно с такими формулами можно ознакомиться в [6, 7, 8, 9, 10]. В то же время следует отметить, что математические модели, базирующиеся на информации о характеристиках уже существующих аппаратов, не могут являться основой для проектирования современных высокоэффективных машин с применением новых технических и технологических решений. Например, с помощью весовых формул трудно учесть влияние силовой схемы крыла или характера распределения воздушной нагрузки по крылу на его массу. В этой связи существует потребность в разработке новых, математически формализуемых методов оценки массы конструкции планера самолёта, использующих в своей основе физические принципы работы конструкции под нагрузкой, описываемые высокоточными численными математическими моделями механики твёрдого тела. В работе [11] представлены научные основы нового метода определения массы конструкций, базирующегося на их силовой работе. Авторы [11] и [12, 13] использует для расчёта теоретической массы силовой конструкции специфический интегральный критерий весового совершенства «силовой фактор» G, характеризующий величину и протяжённость действия усилий в конструкции: , (1) где - эквивалентные напряжения в конструкции, Па, V - объём конструкции, м3 Этот критерий хорошо подходит для поиска и генерации новых силовых схем как с использованием аналитических, так и численных математических моделей механики твёрдого тела. Важнейшим преимуществом такого подхода является возможность оценки весового совершенства именно самой силовой схемы (вне зависимости от её сложности, а также материала конструкции). Критерий «силовой фактор» одновременно учитывает величину и протяжённость действия нагрузки в конструкции, а потому позволяет определять минимально потребный (теоретический) объём и, соответственно, массу равнопрочной полнонапряжённой силовой конструкции, необходимой для восприятия заданного спектра нагрузок с заданным уровнем напряжений. Для проведения весовых оценок нет необходимости подбора сечений силовых элементов: важна лишь силовая схема. В [11, 12] введено также понятие обезразмеренного варианта критерия G - коэффициента силового фактора СК, позволяющего производить сравнение двух или более конкурирующих силовых схем вне зависимости от их абсолютных размеров и величин действующих на них нагрузок. Величина критерия G пропорциональна полной массе исследуемой конструкции mк [11]. В работах [12, 13, 14] приведены примеры эффективного использования критериев G и СК для оценки весовой эффективности силовых схем крыльев и расчёта их массы на основе конечно-элементного моделирования. При этом отдельной задачей является точное определение внешних силовых факторов, в том числе воздушной нагрузки и характера её распределения по площади крыла. Кроме того, во время выбора геометрических параметров крыла самолёта на ранних стадиях проектирования задача параметрической оптимизации его силовой схемы сильно усложняет задачу параметрической оптимизации внешней геометрии крыла за счёт увеличения количества переменных параметров, описывающих его внутреннюю конструкцию. В этой связи целесообразно осуществлять выбор количества и расположения силовых элементов с использованием быстрых и экономных в вычислительном плане аналитических методов оптимизации. Таким образом, целью работы является разработка методики оценки весового совершенства конструкции крыла транспортного самолёта с использованием критерия «коэффициент силового фактора» и малозатратных численных математических моделей механики сплошных сред как эффективного математического инструмента расчёта значения целевой функции при дальнейшей оптимизации параметров геометрического облика крыла и самолёта в целом. 1. УЧЁТ ВЛИЯНИЯ ГЕОМЕТРИЧЕСКИХ ПАРАМЕТРОВ КРЫЛА И СИЛОВОЙ СХЕМЫ НА ЕГО ВЕСОВЫЕ ХАРАКТЕРИСТИКИ На рисунке 1 представлена общая схема функционирования методики, которая предлагается для определения величины критерия СК конструкции крыла с учётом влияния его силовой схемы и внешних геометрических параметров. Предполагается, что в качестве исходных данных выступают следующие величины: - геометрические: площадь крыла, удлинение, сужение, стреловидность, относительная толщина профилей, геометрическая крутка, прочие; - полётные: высота полета и соответствующие ей параметры атмосферы, крейсерская и максимальная скорости полета, максимальная взлетная масса, максимальная эксплуатационная перегрузка. Предусмотрена возможность перевода любой из перечисленных величин из разряда постоянных в разряд проектных переменных, подлежащих варьированию, с целью нахождения их сочетания, сводящего целевую функцию к минимуму. После ввода исходных данных выполняется выбор силовой схемы крыла с определением количества и положения лонжеронов в крыле на основе многокритериальной оптимизации, далее проводится математическое моделирование обтекания крыла на заданном режиме полёта для точного определения величины и характера распределения воздушной нагрузки по площади крыла. Затем выполняется расчёт напряжённого состояния конечно-элементной модели крыла под полученные воздушные нагрузки для расчёта критерия весового совершенства G и его безразмерного коэффициента СК, который выбран в данной работе в качестве эквивалента массы силовой конструкции крыла. Критерий СК, как и G, пропорционален полной массе конструкции [11]. 2. МАТЕМАТИЧЕСКИЕ МОДЕЛИ АЭРОДИНАМИКИ И МЕХАНИКИ ТВЁРДОГО ТЕЛА В качестве математических моделей аэродинамики используется метод дискретных вихрей (МДВ) [15] в линейной стационарной постановке. Данный метод в выбранной постановке использует стационарное несжимаемое потенциальное течение, что, однако не позволяет производить учёт сил вязкости, но даёт возможность с высокой точностью вычислять распределение давления по площади крыла при минимальных вычислительных и временных затратах. Учёт силы вязкости и поправка на сжимаемость среды выполняются с помощью формул инженерной методики расчёта аэродинамических характеристик. Предполагается, что методика рассчитана на самолёты транспортной категории, которые оснащаются крыльями большого удлинения, обладающими в свою очередь линейной зависимостью коэффициента подъёмной силы от угла атаки в диапазоне углов атаки с безотрывным обтеканием. Математическое моделирование по МДВ предполагает вычисление циркуляций в контрольных точках вихревой системы тонкой несущей поверхности (рис. 2 а), геометрически совпадающей со срединной поверхностью крыла, на основе численного решения системы дифференциальных уравнений, отражающих закон Био-Савара-Лапласа: [15, 16]. A = Γ⋅Vn , (2) где A - матрица аэродинамического влияния, 1/м, Γ - матрица-столбец циркуляции вектора скорости, м2/с; Vn - матрица-столбец проекций скорости набегающего потока на нормаль к вихревым рамкам в контрольных точках, м/с. Матрица аэродинамического влияния отражает геометрические характеристики крыла и может быть полностью определена исходными данными. Система уравнений разрешается относительно циркуляций с учётом граничных условий непротекания [15, 16]: , (3) где - проекция скорости набегающего потока на i-ю нормаль; м/с, - вектор набегающего потока, м/с; - нормаль к i-й вихревой рамке; i - номер вихревой рамки. На основе полученного распределения циркуляций рассчитывается элемент нормальной силы в каждой контрольной точке с использованием теоремы Жуковского [15]. На рисунке 2 а представлен общий вид вихревой модели для реализации МДВ. Перед выполнением математического моделирования обтекания крыла выполняется уравновешивание крыла в вертикальной плоскости по соотношению , (4) где - коэффициент подъёмной силы, потребный для обеспечения равенства подъёмной силы весу самолёта с учётом перегрузки, - максимальная эксплуатационная перегрузка, р0 - удельная нагрузка на крыло, qmax - максимальный скоростной напор. При рассмотрении расчётного случая по минимальной эксплуатационной перегрузке в соотношении (4) множитель меняется на min. Затем, путём расчёта по МДВ, выполняется поиск угла атаки a, на котором достигается , потребный для уравновешивания самолёта. Расчёт распределения воздушной нагрузки по крылу выполняется по МДВ для вычисленного угла атаки a. Расчёт напряжённого состояния конструкции крыла и вычисление критерия G, выполняется методом конечно-элементного моделирования в линейной стационарной постановке, для чего автоматически выполняется построение расчётной сетки из оболочечных конечных элементов, воспринимающих растяжение/сжатие, сдвиг и изгиб. Расчёт силового фактора (1) здесь выполняется численно с конечно-элементной модели (КЭМ) (рис. 2 б) по соотношению: , (5) где i - номер элемента модели. Полученный размерный критерий G обезразмеривается до коэффициента СК путём нормирования на характерную нагрузку (вес самолёта m0g) и характерный линейный размер (корень квадратный из площади крыла S): . (6) При рассмотрении нескольких расчётных случаев величина критерия G вычисляется по напряжённому состоянию, полученному из огибающей всех расчётных случаев. 3. ПЕРЕДАЧА НАГРУЗОК ИЗ ВИХРЕВОЙ МОДЕЛИ В КЭМ Отдельной задачей является обеспечение корректного приложения к КЭМ воздушных нагрузок, которые получены аэродинамическим расчётом: размер и взаимное расположение вихревых рамок и конечных элементов в КЭМ могут не совпадать. Кроме того, КЭМ включает в себя нервюры и лонжероны, в совокупности представляющие собой систему жёстких в своей плоскости элементов - стенок, работающих на сдвиг. Однако, расположение некоторых контрольных точек вихревой модели может не совпадать с осью того или иного силового элемента на виде в плане. В таком случае нагрузка, полученная в такой контрольной точке, будет приходиться на нежёсткую из своей плоскости обшивку, что может привести возникновению локальных изгибающих моментов в ней, существенному повышению напряжений и значения силового фактора в этих местах и невозможности адекватно оценить массу конструкции крыла. Для преодоления этих проблем, в работе предложено соблюдать размеры ячеек двух моделей кратными друг другу, например, так, чтобы в одну вихревую рамку на виде в плане вписывались два или четыре элемента КЭМ при совпадении общих узлов. Сила, полученная в контрольной точке вихревой рамки, переносится вдоль линии действия напрямую в соответствующий узел КЭМ, расположенный на верхней поверхности крыла. Кроме того, при построении вихревой модели и КЭМ предварительно высчитывается количество разбиений несущей поверхности на вихревые рамки таким образом, чтобы контрольные точки совпадали с осями силовых элементов на виде в плане и приложение нагрузки приходилось непосредственно на них. 4. ПОИСК ОПТИМАЛЬНОГО РАСПОЛОЖЕНИЯ ЛОНЖЕРОНОВ В КРЫЛЕ ПО ХОРДЕ Задача выбора размещения лонжеронов в крыле по хорде формулируется в терминах нелинейного математического программирования с двумя критериями. В качестве критериев выступают строительная высота лонжеронов с одной стороны и межлонжеронная площадь поперечного сечения крыла - с другой. Строительная высота лонжерона зависит от его положения по хорде, определяется геометрией профиля крыла и подлежит максимизации с целью повышения его изгибной жёсткости. Площадь поперечного сечения крыла, заключённая между лонжеронами, представляет собой замкнутый контур, воспринимающий действие крутящего момента. Она подлежит максимизации для повышения крутильной жёсткости. Проектные переменные - координаты положения лонжеронов по хорде (рис. 3). Рис. 3. Интервалы варьирования проектными переменными При этом полагается, что лонжероны располагаются по процентным линиям крыла, которое обладает прямолинейными передней и задней кромками на виде в плане. В случае рассмотрения двухлонжеронного крыла задача оптимизации формулируется в виде: , (7) где - площадь межлонжеронной части поперечного сечения крыла; - строительная высота 1-го лонжерона; - строительная высота 2-го лонжерона; gв и gН - функции верхней и нижней дужек профиля крыла в сечении; - относительные координаты положения 1-го и 2-го лонжеронов соответственно, отмеренные от передней кромки. Координаты соотнесены к длине хорды данного сечения крыла. Расположение в крыле трёх и более лонжеронов уже необходимо выполнять с учётом распределения крутильных жесткостей между замкнутыми контурами с целью учёта и управления положением оси жёсткости в конструкции. Решение поставленной задачи оптимизации предложено выполнять с использованием многоцелевого эволюционного алгоритма NSGA-II [17], относящегося к группе генетических алгоритмов. 5. РЕАЛИЗАЦИЯ АЛГОРИТМА Изложенная методика алгоритмизирована и функционирует полностью в автоматическом режиме в виде программного обеспечения для ЭВМ на межплатформенной связке двух программных продуктов MATLAB и Ansys Mechaniсal APDL. При этом модуль расчёта аэродинамических характеристик по МДВ, модуль оптимизации положения лонжеронов в крыле, модуль автоматической генерации расчётных сеток для МДВ и КЭМ, а также общая управляющая программа реализованы на языке MATLAB. Расчёт напряжённого состояния КЭМ и вычисление критерия «силовой фактор» - на языке APDL. Блок-схема программы представлена на рисунке 4. 6. РЕШЕНИЕ ДЕМОНСТРАЦИОННОЙ ЗАДАЧИ С целью апробации алгоритма и программного обеспечения выполнено решение прикладной задачи по оценке влияния геометрических параметров крыла гипотетического тяжёлого транспортного самолёта на коэффициент силового фактора его конструкции. Исходные данные: взлётная масса m0=350 т; высота крейсерского полёта Н=10000 м; крейсерская скорость V = 850 км/ч; максимальная эксплуатационная перегрузка = 2,5; минимальная эксплуатационная перегрузка min= -1,5 удельная нагрузка на крыло р0 = 650 ДаН/м2. Геометрические параметры, принятые переменными: стреловидность c, сужение h; постоянные параметры: относительная толщина корневого профиля крыла c0=0,12, концевого - ск=0,08; удлинение l = 8,3; геометрическая крутка j = 0°. На рисунках 5 а, б представлено множество Парето для сочетаний проектных переменных и значений целевых функций при них для профиля крыла с заданными выше геометрическими характеристиками. Из множества сочетаний для двухлонжеронного крыла выбраны относительные координаты положения лонжеронов = 0,23 и = 0,67. Данное сочетание проектных переменных обеспечивает оптимальность по Парето [18]. На рисунках 6 а, б представлены распределения относительной циркуляции по размаху крыла при различных сужениях h и стреловидностях c, полученные с помощью математического моделирования обтекания крыла методом МДВ. На рисунке 7 показана зависимость коэффициента силового фактора СК от этих двух переменных параметров - сужения h и стреловидности c, полученная на основе конечно-элементного моделирования под соответствующие распределения циркуляций. Зависимости, представленные на рисунках 6 - 7, имеют логичный и ожидаемый качественный и количественный характер. Увеличение стреловидности повышает массу крыла из-за увеличения циркуляции в концевых сечениях крыла и уменьшения - в корневых (рис. 6 б). Увеличение сужения - снижает массу за счёт перераспределения воздушной нагрузки ближе к корневым сечениям и снижения изгибающего момента (рис. 6 а). Значения коэффициента силового фактора при соответствующих параметрах количественно и качественно совпадают с результатами, полученными в работах [19, 20], которые посвящены исследованию возможностей применения этого критерия для расчёта массы конструкции крыла транспортных самолётов. Для получения одной точки на графике (рис. 7) по двум расчётным случаям ( и ) потребовалось 25 секунд машинного времени на одном ядре в два потока на процессоре Core i5 с частотой 2,3 ГГц. ЗАКЛЮЧЕНИЕ На основе совместного использования в предложенном алгоритме численных математических моделей высокого уровня из областей механики газа и механики твёрдого тела, а также их реализации в виде программы на ЭВМ получена и апробирована методика определения весовых характеристик крыла большого удлинения. При этом на примере показано, что методика чувствительна к изменению характера распределения воздушной нагрузки по крылу, продиктованной сочетаниями его геометрических параметров, и, следовательно, позволяет учитывать влияние этих параметров на массу конструкции крыла. Это обстоятельство при правильном использовании методики несомненно позволяет определять массу конструкции крыла с более высокой точностью, чем с использованием весовых формул, что оказывает непосредственное влияние на повышение точности определения интегральных показателей эффективности самолёта в целом с сохранением малых временных и вычислительных затрат.

About the authors

O. U Espinosa Barcenas

Samara National Research University

Email: oscar.espinosa.barcenas@gmail.com
Samara, Russia

O. E Lukyanov

Samara National Research University

Email: lukyanovoe@mail.ru
Samara, Russia

References

Supplementary files