Прикладная математика и механика, 2023, T. 87, № 4, стр. 670-683

1. Введение. Прямые численные методы расчета аэродинамики летательных аппаратов (ЛА) требуют существенных компьютерных ресурсов, и даже на современных суперкомпьютерах с алгоритмами распараллеливания для таких расчетов требуются недели и даже месяцы непрерывных вычислений. В связи с этим, начиная с 60-х годов прошлого столетия (когда подобные вычисления были невозможны) широкое применение, в рамках модели идеальной жидкости, получили подходы, основанные на моделировании ЛА в виде вихревых структур, распределенных по поверхности ЛА [1, 2]. С одной стороны, это дает ясную физическую трактовку аэродинамических характеристик, – начиная с простейшей теории несущей линии Прандтля [3]. С другой стороны, это дает надежные численные методы, поскольку такой подход переводит задачу расчета из трехмерной области на двумерные граничные поверхности, что существенно снижает размерность дискретных сеток. В результате сформировался класс методов, которые в российской литературе получили название “методы дискретных вихрей” (методы ДВ) [4–8], а в зарубежной за ними закрепилось название “панельные методы” [9].

Анализ этих методов показывает, что фактически они представляют собой реализацию специальных численных алгоритмов для решения некоторых континуальных одномерных и двумерных интегральных уравнений. Как правило, интегральные уравнения являются сингулярными и/или гиперсингулярными. В указанных выше работах обосновывается корректность перехода от непрерывной к дискретной трактовке в случае одномерных сингулярных ядер, а для одномерных гиперсингулярных ядер обоснование численных алгоритмов основано на их связи с сингулярными. Заметим, что прямое обоснование дискретизации для гиперсингулярных уравнений, без их связи с сингулярными, дано в [10].

Возникающие интегральные уравнения, как одномерные, так и двумерные, решаются для широкого класса геометрий поверхностей, расположенных в потоках жидкости или газа, что после дискретизации приводит к решению систем линейных алгебраических уравнений (СЛАУ) большой размерности, то есть к вычислительной задаче, которая для современных компьютеров не представляет серьезных трудностей.

Исторически развитие численных методов отодвинуло на задний план аналитические подходы как менее точные, поскольку известные аналитические решения основаны, как правило, на некоторых физических допущениях или приближениях. Одним из ярких примеров является классическая линейная теория тонкой несущей поверхности, которая для крыльев большого размаха сводится к одномерному интегро-дифференциальному уравнению Прандтля. Это уравнение самим Прандтлем было получено из физических гипотез, но оно может быть выведено также как теория крыла большого удлинения из двумерного интегрального уравнения линейной теории тонкого крыла [11]. Решение одномерного уравнения Прандтля методом Глауэрта разложением по угловой переменной сводит его к бесконечной СЛАУ [11]. В частном случае незакрученного крыла эллиптической формы в плане и постоянного угла атаки в каждом сечении решение выписывается в явном виде и дает элементарные выражения для аэродинамических характеристик. Считается, что простое приближение для коэффициента подъемной силы ${{c}_{P}} = 2\pi {{\alpha }_{0}}\lambda {\text{/}}(\lambda + 2)$, соответствующее решению Глауэрта для эллиптического крыла, обладает хорошей точностью также для крыльев других форм, при $\lambda \geqslant 4$ [11]. Здесь $\lambda $ – удлинение крыла, ${{\alpha }_{0}}$ – угол атаки. При этом в литературе сложно найти явные решения для крыльев, отличных от эллиптических. Между тем, метод ДВ в его исходной трактовке наиболее естественно приспособлен для крыльев прямоугольной формы в плане, так как именно в этом случае простейшая дискретная сетка по каждой из двух декартовых координат равномерно покрывает всю площадь крыла. Равномерные сетки можно распространить также и на крылья более сложной формы (со скольжением), однако это требует дополнительной трактовки. Авторам неизвестны работы, в которых выражения аэродинамических характеристик для такой естественной канонической геометрии как удлиненное крыло, прямоугольное в плане, были бы представлены в явном виде, основываясь на строгом математическом обосновании. Заполнение этого пробела является одной из целей данной работы. Заодно дается сравнение полученных явных выражений с методом ДВ, а также сравнение с результатами проведенных натурных экспериментов.

2. Асимптотический метод для прямоугольного в плане крыла большого удлинения. Основы метода заложены в [12], где исследуются гармонические по времени колебания прямоугольного в плане крыла, в линейном приближении малых возмущений на фоне набегающего потока. Заметим, что в динамическом случае авторам [12] удалось получить явные аналитические результаты лишь для предельно больших удлинений крыла, в случае так называемого “вырожденного” решения, фактически соответствующего приближению двумерной задачи в каждом сечении крыла вдоль размаха. В рассматриваемом более простом случае стационарного обтекания ниже в рамках трехмерной модели выводятся явные представления для основных аэродинамических характеристик, справедливые в широком диапазоне относительных удлинений крыла.

Пусть тонкое слабоизогнутое прямоугольное в плане крыло размером $( - l,l) \times ( - a,a)$ расположено почти параллельно горизонтальной плоскости xy под малым углом атаки ${{\alpha }_{0}}$ к однородному равномерному потоку, набегающему со скоростью V₀ параллельно оси x, как показано на рис. 1. В общем случае ${{\alpha }_{0}} = - \partial f{\text{/}}\partial x$, где $z = f(x,y)$ – форма поверхности крыла. Далее для простоты в основном рассматривается случай незакрученного крыла с углом атаки, постоянным вдоль оси y.

В рамках линейной теории малых возмущений основное интегральное уравнение метода дискретных вихрей в безразмерной форме имеет следующий вид [4]:

Здесь $\lambda = l{\text{/}}a$ – относительное удлинение крыла, а функция $g(x,y) = \Delta p{\text{/}}(\rho {{V}_{0}})$ связана со скачком давления на крыле $\Delta p = {{p}^{ - }} - {{p}^{ + }}$ при переходе с ее нижней на верхнюю лицевую поверхность, то есть – с локальной подъемной силой, распределенной по поверхности крыла. Заметим, что все величины размерности расстояния приведены к полуширине хорды a.

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

Для асимптотического анализа при $\lambda \to \infty $, представим ядро интегрального уравнения (1.1) в следующем виде:

где интеграл по переменной α понимается как сингулярный интеграл типа Коши.

Будем искать решение интегрального уравнения (2.1) с ядром (2.2) в виде следующего разложения по полиномам Якоби:

которые ортогональны на отрезке (–1, 1) с весом, определяемым качественным поведением решения на передней и задней кромках крыла вдоль хорды. Подстановка соотношений (2.2) и (2.3) в уравнение (2.1) позволяет проинтегрировать по переменной ξ, с использованием табличного интеграла, выражаемого через функции Бесселя [13]:

После этого умножим обе части уравнение (2.1) на функции $\sqrt {(1 + x){\text{/}}(1 - x)} P_{n}^{{\left( { - \frac{1}{2}, + \frac{1}{2}} \right)}}(x)$, $(n = 0,1,...)$ и проинтегрируем по переменной x на интервале (–1, 1). Тогда с использованием табличных интегралов [13]

двумерное уравнение (2.1) преобразуется к бесконечной системе одномерных интегральных уравнений следующего вида: где ${{\delta }_{{nk}}}$ – символ Кронекера.

Как и в [12], в данной стационарной задаче можно показать, что при $\lambda \to \infty $ внедиагональные члены этой системы (при $k \ne n$) асимптотически стремятся к нулю, причем тем быстрее, чем дальше элементы бесконечной матрицы отстоят от главной диагонали. Более точная оценка получается после перехода на отрезок (–1, 1) по переменным $y,\eta $ заменой переменных $\tilde {y} = y{\text{/}}\lambda $, $\tilde {\eta } = \eta {\text{/}}\lambda $, $\tilde {\beta } = \lambda \beta $. После этого явная асимптотическая оценка ядер становится такой:

Следовательно, главный член асимптотики можно получить, сохраняя лишь диагональные члены (при $k = n$). Таким образом, бесконечная система одномерных интегральных уравнений распадается на множество одномерных интегральных уравнений, которые можно решать независимо друг от друга. При этом в общем случае произвольной формы крыла $f(x,y)$ получаем бесконечное число независимых одномерных интегральных уравнений. В частном случае постоянного угла атаки $\partial f(x,y){\text{/}}\partial x$ = = $ - {{\alpha }_{0}} \equiv \operatorname{const} $ правая часть ${{F}_{n}}(y)$ для всех $n$, кроме $n = 0$, равна нулю: ${{F}_{0}}(y) = {{V}_{0}}{{\alpha }_{0}}$, ${{F}_{n}}(y) \equiv 0$ $(n \geqslant 1)$. В результате, все неизвестные функции ${{g}_{k}}(y)$, кроме ${{g}_{0}}(y)$, асимптотически равны нулю, а для ${{g}_{0}}(y)$ получаем следующее интегральное уравнение:

При этом из (2.3) следует, что решение уравнения (2.1) принимает вид

Асимптотика решения уравнения (2.7) при большом удлинении λ основана на реализации идеи построения погранслойных решений. В центральной зоне крыла, вдали от боковых кромок, концы интервала интегрирования по переменной η можно удалить на бесконечность. Построение такого “вырожденного” решения основано на решении уравнения в свертках и легко приводит к явному выражению применением преобразования Фурье (ПФ) по у:

Здесь функция $G_{0}^{c}(\beta )$ обозначает ПФ от функции $g_{0}^{c}(y)$, а $\delta (\beta )$ – дельта-функция Дирака. Очевидно, что вырожденное решение справедливо в “центральной зоне крыла” (вдали от боковых кромок, отсюда индекс “c” в соответствующих функциях) и соответствует плоской задаче. При этом, например, подъемная сила P и коэффициент подъемной силы ${{c}_{P}}$ равны соответственно

что действительно соответствует двумерной теории [3, 11].

Структура решения уравнения (2.7) в окрестности левой кромки $y = - \lambda $ может быть получена удалением правой кромки $y = \lambda $ на бесконечность. В результате приходим к уравнению Винера–Хопфа (В–Х)

записанному в координатах, отсчитываемых от левой кромки $y = - \lambda $. Такой подход к построению погранслойного решения был предложен в контактных задачах теории упругости [14]. Родственные подходы развивались в волновых задачах [15, 16], где похожий метод называется “метод краевых волн”.

Метод В–Х подробно описан в литературе [14–16]. В применении к уравнению (2.10) основная сложность состоит в факторизации функции ${{L}_{0}}(\beta )$, которую решаем методом приближенной факторизации. С этой целью приблизим эту функцию на вещественной оси так, чтобы учесть ее поведение в нуле и на бесконечности и обеспечить простую факторизацию. Поскольку

то простейшее приближение обеспечивает точность с относительной погрешностью не более 17% во всем диапазоне $0 \leqslant \left| \beta \right| < \infty $ и гарантирует простую факторизацию при применении метода В–Х. Технические детали решения уравнения (2.10) с символьной функцией ${{L}_{0}}(\beta )$ в виде (2.12) довольно просты, см. например [15, 16]. В результате для ПФ решения уравнения Винера–Хопфа (Wiener–Hopf) получаем где переход к переменной $p = - i\beta $ позволяет от обратного ПФ от комплексной функции перейти к обратному преобразованию Лапласа от вещественной функции, что с использованием таблиц [17] приводит к следующему решению выражающемуся через функцию вероятности ${\text{Erf}}\left( t \right)$.

При этом, согласно методу малого параметра [14–16], главный член асимптотики решения уравнения (2.7) при больших удлинениях $\lambda $ может быть представлен в виде:

что с учетом (2.8) дает полное представление решения.

К сожалению, погрешность аппроксимации (2.12) 17% не позволяет достичь желаемой точности в сравнении с расчетом методом ДВ. В связи с этим подбирается более точное приближение:

где положительные константы $\gamma ,\mu ,\nu $ выбраны из условия наилучшего приближения. Расчеты показывают, что при подобранных значениях констант, указанных в (2.16), обеспечивается точность с относительной погрешностью, не превышающей 3.5% на всем интервале $0 \leqslant \left| \beta \right| < \infty $. При этом решение уравнения В–Х в образах Лапласа ($p = - i\beta $) принимает вид: что после обращения преобразования Лапласа приводит к следующему решению [17]:

Заметим, что, поскольку $\mu < \gamma $, то вторая функция вероятности в (2.18) имеет мнимый аргумент.

3. Вычисление аэродинамических характеристик крыла. Любопытно вычислить подъемную силу на основе представления (2.18):

Тогда для безразмерного коэффициента подъемной силы получаем:

Здесь учтены табличные интегралы [18]:

Собирая подобные члены в (3.2), приходим к следующему выражению:

Здесь учтено, что при $\sqrt {\mu - \gamma } = i\sqrt {\gamma - \mu } = i\sigma $, $(\sigma > 0)$ имеет место следующее соотношение:

где ${\text{Erfi}}\left( t \right) = \left( {2{\text{/}}\sqrt \pi } \right)\int_0^t {\exp \left( {{{\tau }^{2}}} \right)d\tau } $ – функция вероятности мнимого аргумента [19, 20].

Для больших удлинений $\lambda $ выражение (3.3) может быть упрощено с использованием следующих асимптотических представлений функций вероятности для большого аргумента:

где величинами 5-го порядка малости пренебрегаем. Тогда формула (3.3) примет следующий вид:

Поскольку при сложении показателей двух последних экспоненциальных функций имеем $1{\text{/}}\mu - 1{\text{/}}\delta = 1{\text{/}}\gamma $, то (3.5) можно переписать в виде:

или, подставляя значения числовых параметров, в окончательном виде:

Заметим, что первый член во второй строке асимптотического разложения (3.5), убывающий как квадратный корень из удлинения, сократился после применения разложений (3.4).

Сравнение с классическим приближением по формуле Глауэрта [11] c_P = $2\pi {{\alpha }_{0}}\lambda {\text{/}}(\lambda + 2)$, а также с решением по методу ДВ приведено в табл. 1. Заметим, что первые два члена асимптотики при большом удлинении в формуле Глауэрта имеют вид: ${{c}_{P}}{\text{/}}{{\alpha }_{0}}$ = = $2\pi (1 - 2{\text{/}}\lambda )$, в то время как в формуле (3.7): ${{c}_{P}}{\text{/}}{{\alpha }_{0}}$ = $2\pi (1 - 2.401{\text{/}}\lambda )$. Таким образом, выход на решение плоской задачи ${{c}_{P}}{\text{/}}{{\alpha }_{0}}$ = $2\pi = 6.283$ с ростом удлинения для прямоугольного крыла происходит медленнее, чем для эллиптического. Это также видно по различию погрешностей этих решений в сравнении с методом ДВ в табл. 1 для существенно больших удлинений.

Заметим, что в литературе имеется много полуэмпирических формул для подъемной силы прямоугольных крыльев, более точных в сравнении с экспериментом, чем формула Глауэрта [11]. Применяются также комбинации полуаналитических и эмпирических подходов (смотри, например, [21]). Однако из известных аналитических приближений лишь теория Прандтля–Глауэрта, отраженная в табл. 1, основана на теоретическом фундаменте с достаточно строгим обоснованием.

Перейдем к вычислению индуктивного сопротивления. Для вычисления этой величины известны различные формулы, выражающиеся через интегралы от произведения циркуляции и вертикальной компоненты вектора скорости потока в плоскости вихревой пелены [1, 22]. Здесь применим альтернативный, но достаточно естественный метод, при котором индуктивное сопротивление может быть вычислено как проекция подъемной силы на направление потока минус подсасывающая сила. В линейном приближении малого угла атаки имеем:

где подъемная сила P легко вычисляется через безразмерный коэффициент (3.7), а выражение для подсасывающей силы T известно в теории тонкого крыла [23]. Детальный анализ явных представлений в формулах (3.1) и (3.8) показывает, что выражение ${{D}_{{{\text{ind}}}}}$ в (3.8) всегда положительно и стремится к нулю при $\lambda \to \infty $, что соответствует выходу на плоскую задачу, как и должно быть по физическому смыслу. В общем случае интеграл в (3.8) вычислялся по квадратурной формуле.

4. Сравнение с экспериментальными данными. Для сравнения с теоретическими результатами были проведены натурные эксперименты по обдуву в аэродинамической трубе. Образец представляет собой прямоугольную дюралюминиевую пластинку толщиной 2 мм габаритами 20 × 2.5 см в плане, с удлинением $\lambda = 20{\text{/}}2.5$ = 8. Обдув осуществлялся в аэродинамической трубе замкнутого типа компании “Денар” (г. Ярославль) с рабочей частью 30 × 30 см в поперечнике и длиной 60 см. Скорость набегающего потока 11 м/с. Подъемная сила P измерялась при помощи аэродинамических весов, скорость набегающего потока – при помощи трубок Пито–Прандтля. Угол атаки ${{\alpha }_{0}} = - \partial f{\text{/}}\partial x$ постоянный как вдоль размаха, так и вдоль хорды в каждом сечении (прямолинейная хорда). В эксперименте ${{\alpha }_{0}}$ изменялся от 2 до 16° с шагом 2°. Стабильность результатов при повторных измерениях – от 10 до 15%.

Оценим режим обтекания в натурном эксперименте – ламинарный он или турбулентный. Число Рейнольдса равно $\operatorname{Re} = 2a{{V}_{0}}{\text{/}}{{\nu }_{{\text{в}}}}$ = (0.025 м) (11 м/с)/($1.5 \times {{10}^{{ - 5}}}$ м²/с) = = 18 333, что, по крайней мере, на порядок меньше критического числа $\operatorname{Re} {\kern 1pt} * = 3.5 \times {{10}^{5}}$ при переходе от ламинарного режима к турбулентному при обтекании тонкой пластинки [24]. Таким образом, в эксперименте – ламинарный режим. Здесь ${{\nu }_{{\text{в}}}}$ – кинематическая вязкость воздуха.

Для аэродинамического качества Q, с целью сравнения теоретических и экспериментальных результатов, добавим к индуктивному сопротивлению (3.8) силу вязкого трения, известную при продольном обтекании тонкой пластинки из теории Прандтля для ламинарного погранслоя как формула Блазиуса [24]:

которую для малого угла атаки считаем не зависящей от угла ${{\alpha }_{0}}$. Сравнение теоретических и экспериментальных результатов для аэродинамического качества Q = P/D отражено в табл. 2.

Заключение. 1. В работе предлагается аналитическая теория прямоугольного крыла в рамках гипотезы малых линейных возмущений. Она является асимптотической по большому удлинению крыла, однако основной результат, представленный формулой (3.7), показывает на сравнении с расчетами по методу ДВ (табл. 1), что ее точность является хорошей также и в области средних удлинений.

2. Считается, что формула Глауэрта для коэффициента подъемной силы ${{с}_{P}} = 2\pi {{\alpha }_{0}}\lambda {\text{/}}(\lambda + 2)$ является достаточно универсальной для крыльев произвольной формы в плане (без стреловидности) и обладает хорошей точностью при $\lambda \geqslant 4$. Табл. 1 позволяет оценить ее точность в применении к прямоугольному крылу.

3. Сравнение теории и эксперимента по аэродинамическому качеству, отраженное в табл. 2, позволяет сделать любопытные выводы. Две колонки с теоретическими результатами, соответственно 3-я и 4-я колонки, отличаются следующим. В первом случае в знаменателе – полное сопротивление D, равное проекции подъемной силы на направление потока: $P{{\alpha }_{0}}$, минус подсасывающая сила T, плюс вязкое сопротивление W. Во втором случае в знаменателе не учитывается подсасывающая сила T, и оказывается, что это гораздо ближе к экспериментальным данным, чем результаты 3-й колонки. С физической точки зрения это вполне закономерно, так как подсасывающая сила возникает при обтекании бесконечно-тонкой пластинки, которую трудно считать таковой при толщине 2 мм и длине хорды 2.5 см. Фактически в эксперименте никакой подсасывающей силы нет, и реальное сопротивление создается подъемной силой, плюс вязкое трение. Сравнение теории из 4-й колонки с экспериментом графически представлено на рис. 2. Как обычно, теория дает несколько завышенные значения качества Q, что на рисунке проявляется практически при всех углах атаки.

4. Сравнение 3-й и 5-й колонок в табл. 2. Известно, что среди всех крыльев одного и того же удлинения наименьшее индуктивное сопротивление имеет крыло эллиптической формы в плане [3, 11, 21]. Из табл. 2 следует, что это свойство, в частности, проявляется в том, что аэродинамическое качество Q у эллиптического крыла при любом угле атаки больше, чем у прямоугольного.

5. Результат, представленный в формуле (2.8), может показаться равносильным тому, чтобы изначально искать решение в виде произведения двух функций, зависящих соответственно от x и от y, т.е. в виде функции с разделяющимися переменными. На самом деле такой подход не эквивалентен предложенному в данной работе по следующим причинам.

Во-первых, если искать асимптотическое решение двумерного уравнения (2.1) в виде (2.8), то непонятно как из двумерного уравнения получить одномерное уравнение для нахождения функции ${{g}_{0}}(y)$. Например, если приравнять левые и правые части в (2.1) при каком-то одном значении переменной x (как простейший вариант метода коллокации), то полученное одномерное уравнение не будет эквивалентно уравнению (2.7); как следствие – не удастся обосновать, что полученное таким образом решение будет давать главный член асимптотики при больших удлинениях крыла.

В связи с отмеченным в предыдущем абзаце заметим, что в [11] при выводе одномерного интегро-дифференциального уравнения Прандтля применяется скалярное умножение (по переменной x) уравнения (2.1) на функцию $\sqrt {(1 + x){\text{/}}(1 - x)} $. При этом не делается никаких асимптотических оценок относительно связи полученного одномерного уравнения с исходным двумерным (2.1). В отличие от этого, в настоящей работе обосновывается, что в рамках предложенного подхода двумерное интегральное уравнение асимптотически распадается на независимые одномерные уравнения, что и позволяет построить решение в виде (2.8) + (2.18).

Работа выполнена в рамках проекта РФФИ, грант № 19-29-06013.

Авторы также признательны Грунтфесту Р.А. за внимание к работе.