Modeling of the unsteady aerodynamic characteristics of the NACA 0015 airfoil from the data of numerical calculations of the flow

全文:

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

Как правило, идентификация параметров математических моделей осуществляется по данным динамических экспериментов [4, 5], реже по летным данным [6, 7]. В первом случае известным ограничением является невозможность обеспечить подобие по числу Рейнольдса во время испытаний в аэродинамических трубах, во втором — ограниченное количество летных данных на больших углах атаки и сложность их первичной обработки. Численные методы лишены указанных недостатков, поэтому их использование для идентификации математических моделей нестационарной аэродинамики представляет значительный интерес.

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

В литературе исследуется рациональность применения существующих методов и моделей в различных задачах [8−13]. В работах [14−16] рассматривались как классические подходы, решение уравнений Рейнольдса в нестационарной постановке с различными моделями турбулентности (URANS), так и гибридные URANS—LES-методы (методы моделирования отсоединенных вихрей (DES, DDES)), методы крупных вихрей (LES). Лучшее согласование достигается при использовании LES-подхода. Однако при этом затраты на проведение таких расчетов самые высокие по сравнению с другими указанными методами.

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

В настоящей работе на примере нестационарного обтекания профиля NACA 0015 при колебаниях по углу атаки с различными частотами и амплитудами исследуется возможность применения результатов численного моделирования для разработки приближенной математической модели нестационарной аэродинамики, пригодной для использования в задачах проектирования роторов ветровых турбин с горизонтальной и вертикальной осью, а также роторов вертолетов. Для этого решаются уравнения Рейнольдса в стационарной и нестационарной постановках с двумя моделями турбулентности: Спаларта–Альмараса (S–A) и k–ω SST. Проводится валидация результатов расчета путем сравнения с данными эксперимента. По этим данными идентифицируется модель нестационарных коэффициентов подъемной силы и продольного момента, действующих на профиль при его неустановившемся движении на больших углах атаки. Модель формулируется в рамках подхода с введением внутренней динамической переменной [3]. Результаты моделирования сравниваются с расчетными и экспериментальными данными.

1. ЧИСЛЕННЫЙ РАСЧЕТ

Выбор эксперимента для валидации численного расчета

Колебания профиля NACA 0015 исследовались в экспериментах [17−19]. Для выбора численной постановки задачи исследования обтекания колеблющегося по углу атаки профиля крыла проведено сравнение с экспериментальными данными работы [18], которые являются тестовым случаем для многих численных работ.

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

Исследовалась модель прямого крыла с профилем NACA 0015 с хордой c = 0.3048 м в сечении, полуразмах модели L = 1.524 м. В работе есть часть экспериментальных исследований, где для обеспечения двумерности течения у боковых сечений крыла устанавливались специальные пластины. Это позволит далее использовать двумерную постановку численной задачи.

В качестве турбулизаторов потока использовалась зубчатая лента на передней кромке крыла. По размаху крыла в нескольких сечениях располагались датчики давления. В данной работе приводится сравнение экспериментальных аэродинамических характеристик в эксперименте в двух сечениях: y/L = 0.263 (10 дифференциальных датчиков) и y/L = 0.500 (20 датчиков давления). Для интегрирования проводилась интерполяция распределения давления в сечении по специальной процедуре. Соответственно, данные о силах и моментах из расчета для сравнения с экспериментом приводятся только для компоненты давления.

Постановка численной задачи

Рассматривается режим обтекания профиля крыла NACA 0015 с числом Маха M = 0.29 и числом Рейнольдса Re = 1.94 · 10⁶ в набегающем потоке. Угол атаки изменяли в диапазоне α = 0−20°. В расчетах для верификации и валидации заданы точные параметры из эксперимента. Геометрическая модель профиля соответствует используемому в эксперименте.

Вокруг профиля была построена структурированная многоблочная расчетная сетка C—H-топологии. Сеточные сгущения сделаны в области по нормали к поверхности профиля (у первой ячейки Y₊ ~ 1), а также в области задней кромки и за профилем в области следа.

Для получения обтекания стационарного профиля решаются двумерные уравнения Рейнольдса для сжимаемого газа в стационарной постановке (2D RANS), для случая с колебаниями — в нестационарной постановке (2D URANS). Для замыкания уравнений используются различные модели турбулентности.

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

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

$\alpha ={\alpha }_0-\Delta \alpha \sin 2\pi ft,$

где α₀ — средний угол атаки, Δα — амплитуда колебаний, f — частота колебаний, t – время.

Для каждого набора параметров рассчитывали несколько периодов колебаний, а для сравнения с экспериментальными данными выбирали период, где течение уже установилось.

Верификация

Для определения размера расчетной сетки, на которой возможно получить результаты, требуется проверить сходимость результатов по сеткам. Для этого рассматриваются различные размеры расчетных сеток (70 000, 250 000, 460 000, 1 100 000 ячеек). Сетки были получены увеличением количества ячеек в двух направлениях.

Проведены исследования обтекания стационарного профиля крыла. В данном блоке расчетов для замыкания уравнений Рейнольдса использовали однопараметрическую модель турбулентности S–A.

Полученные результаты показывают, что, начиная с сетки размером 250 000 ячеек, решение не меняется. Соответственно, для проведения дальнейших расчетов будет использоваться сетка такого размера.

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

Для режима с M = 0.29 и α = 20° были рассмотрены 2 шага по времени Δt = 2·10^–6, 10^–5 и несколько различных точностей сходимости Δ = 5 · 10^–6, 10^–5, 5 ·10^–5, 5 ·10^–4. В результате указанных расчетов было определено, что можно проводить исследования с точностью сходимости Δ = 5 ·10^–4 и шагом по времени Δt = 10^–5, так как полученные кривые не различаются в пределах точности определяемых параметров.

Валидация

Для сравнения с экспериментальными данными рассматривались две модели турбулентности — S–A- и k–ω-SST. На рис. 1 приведены аэродинамические характеристики для указанных моделей и экспериментальные данные для стационарного профиля в двух сечениях крыла.

 

Рис. 1. Сравнение зависимости коэффициента нормальной силы от угла атаки C_y(α)-профиля для экспериментальных и численных результатов с различными моделями турбулентности для стационарного профиля крыла: 1 — эксперимент, сечение y/L = 0.5, 2 — эксперимент, сечение y/L = 0.263, 3 — модель турбулентности S–A, 4 — модель турбулентности k–ω-SST; M = 0.29.

 

Результаты расчетов с моделью k–ω-SST (рис. 1) лучше согласуются с экспериментальными данными для сечения y/L = 0.500 по C_y, чем модель S–A, в особенности в отрывной области при α > 16°. На линейном участке и в области максимального коэффициента нормальной силы сдвиг, по-видимому, обусловлен влиянием трехмерности течения. Сдвиг по углу атаки для модели SST меньше, чем для модели S–A, и составляет Δα ≈ 1.5°.

На рис. 2 представлены кривые для случая с колебаниями профиля по углу атаки при α₀ = 10.88°, Δα = 4.22°, f = 10.1 Гц.

 

Рис. 2. Сравнение зависимости коэффициента нормальной силы C_y(α) от угла атаки профиля для экспериментальных и численных результатов с различными моделями турбулентности для колеблющегося профиля α₀ = 10.88°, Δα = 4.22°, f = 10.1 Гц: 1 — эксперимент, сечение y/L = 0.5, 2 — эксперимент, сечение y/L = 0.263, 3 — модель турбулентности S–A, 4 — модель турбулентности k–ω-SST; M = 0.29.

 

По форме зависимости C_y(α) и значениям модель k–ω-SST лучше описывает экспериментальную кривую. Для прямого хода по углам атаки, начиная с α = 10°, m_z(α) лучше моделируется с помощью модели S–A, но начальный участок прямого хода и конечный обратного лучше при выборе k–ω-SST.

Для проведения расчетов, таким образом, далее будет использоваться модель турбулетности k–ω-SST.

Результаты расчетов

Основываясь на выбранных параметрах задачи, проведены параметрические расчеты колебаний профиля: вычисления проводились на сетке размером 250 000 ячеек, для замыкания уравнений использовалась модель турбулетности k–ω-SST, а шаг по времени выбирался исходя из сохранения количества точек на период колебаний.

Параметры набегающего потока: число Маха M = 0.1428, число Рейнольдса Re = 0.97·10⁶. Проведено два блока расчетов с частотами f = 0.5, 1, 2 Гц: расчеты колебаний профиля со средними углами атаки α₀ = 0−16° с шагом Δα₀ = 4° и амплитудой Δα = 4° (рис. 3) и расчеты с большей амплитудой Δα = 8° на средних углах атаки α₀ = 4, 11, 15° (рис. 4). Дополнительно проведен расчет с α₀ = 11°, Δα = 9°. Результаты данных расчетов далее будут использованы для построения математической модели.

 

Рис. 3. Сравнение зависимости коэффициента нормальной силы от угла атаки C_y(α) и коэффициента момента тангажа m_z(α) профиля для численных результатов: 1 — обтекание стационарного профиля, колебания c амплитудой Δα = 4° частотой f = 2 Гц и различными средними углами атаки α₀: 2 — α₀ = 0, 3 — α₀ = 4°, 4 — α₀ = 8°, 5 — α₀ = 12°, 6 — α₀ = 16°; M = 0.1428.

 

Рис. 4. Сравнение зависимости коэффициента нормальной силы от угла атаки C_y(α) и коэффициента момента тангажа m_z(α) профиля для численных результатов при a, б –α₀ = 4; в, г — α₀ = 11°: 1 — обтекание стационарного профиля, колебания с амплитудой Δα = 8 и различными частотами: 2 — f = 0.5 Гц, 3 — f = 1 Гц, 4 — f = 2 Гц; M = 0.1428.

 

По рис. 3а видно, что при α < 12° ширина гистерезисной петли практически одинакова для всех средних углов атаки. При α ≥ 12° увеличение среднего угла атаки приводит к увеличению ширины гистерезисной петли.

При увеличении частоты колебаний (рис. 4а, 4в) ширина гистерезисной петли также растет. Значение C_ymax в случае динамического гистерезиса выше, чем в случае обтекания статического профиля, и достигается позже по углу атаки. Для m_z зависимость от частоты аналогичная (рис. 4б). На большем среднем угле атаки при увеличении частоты колебаний ширина кривой гистерезиса не только увеличивается, но и изменяется вид кривой в области α₀ ≈ 18° (рис. 4г).

2. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

Структура математической модели аэродинамики

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

$C=f\left(\text{ξ\hspace{0.17em}},x\right),\frac{dx}{dt}=g\left(\text{ξ\hspace{0.17em}},x\right).$

Здесь C = (C_y, m_z), $\xi$ = (α, ω_z) — вектор кинематических параметров движения, x — вектор внутренних переменных, определяющих состояние отрывного обтекания.

Для профиля на больших углах атаки внутренней переменной можно считать, например, положение точки отрыва потока x_s вдоль хорды профиля. В линейной постановке в качестве функции f может быть выбрано решение задачи Кирхгофа об отрывном обтекании профиля идеальной несжимаемой жидкостью [20, 21].

В рамках этой задачи считается, что отрыв потока происходит в некоторой точке верхней поверхности и с задней кромки, линии срыва потока, распространяясь вниз по течению, остаются вблизи оси Ox до бесконечности. Численное сравнение решения этой задачи с решением нелинейной задачи Чаплыгина–Лаврентьева с обтеканием пластинки со срывом струи с верхней поверхности [22], в рамках которой не ставится ограничивающие предположение о форме срывной зоны, показывает, что при углах атаки α < 20° качественные различия отсутствуют, и решение Кирхгофа для коэффициента нормальной силы

$C_y\left(\alpha ,x_s\right)=\frac{\pi }{2}\sin \alpha {\left(1+\sqrt{x_s}\right)}^2$ (2.1)

можно считать приближенно верным. Учет конечности зоны отрыва потока [23] также не приводит к существенным качественным изменениям.

Таким образом, формулу Кирхгофа (2.1) можно использовать для моделирования нестационарного поведения коэффициента C_y при неустановившемся движении профиля. При тех же предположениях можно вывести и формулу для коэффициента момента тангажа, но оценки показывают, что она менее адекватна. На рис. 5 представлены красными сплошными линиями с маркерами результаты статических расчетов методами CFD-зависимостей C_y(α) (верхний график) и m_z(α) (нижний график).

 

Рис. 5. Зависимости С_y, m_z и x_s от α. 1 — Решение Кирхгофа, 2 — статические CFD-расчеты (маркеры), 3 — аппроксимация Кирхгофа, 4 — аппроксимация при безотрывном обтекании, 5 — аппроксимация при полностью отрывном обтекании.

 

Если рассматривать формулу (2.1) как уравнение для нахождения x_s(α) по известной зависимости C_y(α), можно получить предполагаемую зависимость положения точки отрыва на верхней поверхности профиля x_s(α). Эта зависимость изображена на среднем графике рис. 5 черной пунктирной линией. При взятии этой стационарной зависимости за основу в предположениях Кирхгофа зависимости для коэффициентов подъемной силы и момента тангажа принимают вид: зеленые сплошные линии на верхнем графике − C_y(α) и соответствующая линия на нижнем графике — m_z(α). Видно, что для C_y(α) аппроксимация весьма удовлетворительна, тогда как для m_z(α) оставляет желать лучшего. Поэтому математическая модель для подъемной силы строилась на основе формулы (2.1) с дополнительным динамическим уравнением, описывающим запаздывание развития отрыва потока для внутреннего переменного x_s(α), а для момента тангажа изменение m_z предлагается определять через изменение подъемной силы и положение фокуса профиля Δx(α):

$m_z(\alpha ,\Delta x)=m_{z_0}+\Delta x\left(\alpha \right)C_y\left(\alpha \right).$ (2.2)

Здесь $m_{z_0}$ — в общем случае ненулевой момент профиля при нулевой подъемной силе. В рассматриваемом случае в связи с симметричностью профиля NACA 0015 $m_{z_0}=0.$

Таким образом, второй внутренней динамической переменной математической модели является Δx. В этом случае предлагаемая статическая аппроксимация m_z(a) совпадает с результатами CFD-расчета удовлетворительно (см. рис. 5, нижний график, черная пунктирная линия).

Учитывая сказанное, математическая модель продольных аэродинамических характеристик профиля может быть записана в виде

$\begin{array}{c}{C}_y\left(t\right)=C_y\left(\alpha ,x_s\right)+C_y^{{\overline{\omega }}_z}{\overline{\omega }}_z,m_z\left(t\right)=C_y\left(t\right)\Delta x\left(t\right)+m_z^{{\overline{\omega }}_z}{\overline{\omega }}_z,\\ {\tau }_1{\dot{x}}_s+x_s=x_{s_0}\left(\alpha -{\tau }_2\overline{\dot{\alpha }}\right),{\tau }_3\Delta \dot{x}+\Delta x+k_1\Delta x^2+k_2\Delta x^3=\Delta x_0\left(\alpha -{\tau }_4\overline{\dot{\alpha }}\right).\end{array}$ (2.3)

Здесь $C_y^{{\overline{\omega }}_z}$ и $m_z^{{\overline{\omega }}_z}$ — вращательные производные, описывающие эффекты демпфирования присоединенного потока; $x_{s_0}$ — положение точки отрыва потока в статических условиях, которое однозначно определяется по численной статической зависимости C_y(α) при помощи формулы (2.1) (см. рис. 5, средний график); τ₂, τ₄ — безразмерные характеристические времена запаздывания отрыва потока, τ₁, τ₃ — безразмерные характеристические времена релаксационных процессов динамической системы, Δx₀ — изменение аэродинамического фокуса в статических условиях [3]. В дифференциальное уравнение для Δx включены также нелинейные слагаемые с коэффициентами k₁, k₂.

Таким образом, сформулирована математическая модель (2.3) в виде двух нелинейных соотношений и двух обыкновенных дифференциальных уравнений, описывающих эффекты запаздывания. Модель для продольных аэродинамических характеристик профиля содержит восемь неизвестных параметров τ₁, τ₂, τ_3, τ_4, k₁, k₂, $C_y^{{\overline{\omega }}_z}$, $m_z^{{\overline{\omega }}_z}$, которые необходимо идентифицировать по результатам численного расчета при вынужденных колебаниях профиля.

Идентификация параметров модели

Идентификация проводилась на части данных численного расчета вынужденных гармонических колебаний, а именно рассматривались колебания около среднего угла атаки α₀ = 11° с амплитудой Δα = 8° и частотой 2 Гц, осредненных на один период путем минимизации функционалов рассогласования численного и смоделированного значений на заданном массиве точек полного периода длиной n

${\Phi }_{C_y}({\tau }_1,{\tau }_2,C_y^{{\overline{\omega }}_z})=\sum _{i=1}^n{\left(C_y^{CFD}\left(t_i\right)-C_y^{\mathrm{mod}}\left(t_i\right)\right)}^2\text{},{\Phi }_{m_z}\left({\tau }_3,{\tau }_4,k_1,k_2,m_z^{{\overline{\omega }}_z}\right)=\sum _{i=1}^n{\left(m_z^{CFD}\left(t_i\right)-m_z^{\mathrm{mod}}\left(t_i\right)\right)}^2.$

Были получены следующие параметры. Для С_y: τ₁ = 2.31, τ₂ = 4.32, $C_y^{{\overline{\omega }}_z}$ = -7.46, для m_z: τ₃ = 0.1, τ₄ = 0.69, $m_z^{{\overline{\omega }}_z}$= -1, k₁ = -6.17, k₂ = -20.34. На рис. 6 на верхних графиках маркерами показаны результаты расчетов, по которым проводилась идентификация, сплошными синими линиями — результаты моделирования. На нижних графиках — соответствующие динамические положения точки отрыва и аэродинамического фокуса. Пунктирными линиями на всех графиках приведены статические характеристики C_y(α, x_s0), m_z(α, x_s0), x_s0, Δx₀. Для коэффициента нормальной силы согласование удовлетворительное, для коэффициента момента тангажа модель занижает значение m_z в диапазоне α = 12−17° при увеличении угла атаки и завышает в диапазоне α = 10−17° при уменьшении α. Зелеными сплошными линиями на графиках показаны C_y(t) и m_z(t) без учета линейных слагаемых $C_y^{{\overline{\omega }}_z}$, $m_z^{{\overline{\omega }}_z}.$ При малых углах атаки, α < 10°, значительная петля динамического гистерезиса в m_z обусловлена большим в относительных величинах линейным демпфированием. Самопересечение петли при α = 16° говорит о смене знака демпфирования. Динамические значения C_y при безотрывном обтекании близки к статическим.

 

Рис. 6. Результаты идентификации параметров математической модели, α₀ = 8°, Δα = 11°, f = 2 Гц. 1 — статика, 2 — без учета линейного слагаемого, 3 — с учетом линейного слагаемого, 4 — результаты CFD расчета, 5 — статическое изменение внутренней переменной, 6 — динамическое изменение внутренней переменной.

 

Результаты моделирования

При помощи разработанной математической модели продольных аэродинамических характеристик, идентифицированной по части результатов CFD, проводилось моделирование нестационарных аэродинамических нагрузок при гармонических колебаниях профиля для всех расчетных случаев. На рис. 7–10 показаны результаты моделирования (линии) для нескольких случаев в сравнении с расчетом (маркеры).

 

Рис. 7. Динамическое изменение C_y, m_z, x_s, Δx при гармонических колебаниях с α₀ = 8°, Δα = 11° и частотами f = 0.5, 1 Гц.

 

Рис. 8. Динамическое изменение C_y, m_z, x_s, Δx при гармонических колебаниях с α₀ = 12°, Δα = 4° и частотами f = 0.5, 1, 2 Гц.

 

Рис. 9. Динамическое изменение C_y, m_z, x_s, Δx при гармонических колебаниях с α₀ = 16°, Δα = 4° и частотами f = 0.5, 1, 2 Гц.

 

Рис. 10. Динамическое изменение C_y, m_z, x_s, Δx при гармонических колебаниях с α₀ = 11°, Δα = 9° и частотами f = 0.5, 1, 2 Гц.

 

Пунктирные линии на всех графиках соответствуют статическим зависимостям нормальной силы, момента тангажа, положения точки отрыва и аэродинамического фокуса от угла атаки. Для C_y согласование результатов моделирования с расчетными данными удовлетворительное для всех частот колебаний как при безотрывном обтекании (см. рис. 8), так и в области α = 14−20°, где отрывное обтекание обусловливает возникновение значительных петель динамического гистерезиса (см. рис. 6, 9). Динамическое положение точки отрыва во всех случаях существенно отличается от своего статического значения. Можно отметить отличие расчетной кривой для частоты f = 2 Гц в области смены направления движения, связанное, по-видимому, с эффектами, не учитываемыми в модели.

Для коэффициента момента тангажа результаты удовлетворительно согласуются для частот f = 0.5, 1 Гц, для f = 2 Гц согласование носит скорее качественный характер, модель занижает значение m_z в области α ≈ 14−18°.

При большой частоте f ≈ 10 Гц тенденции аналогичные. Результаты расчета и моделирования m_z близки при колебаниях около α₀ ≈ 11° с амплитудой 4° (рис. 11), для α₀ ≈ 15° и той же амплитуды расчет демонстрирует более широкие гистерезисные петли.

 

Рис. 11. Динамическое изменение C_y, m_z при гармонических колебаниях: а) α₀ = 10.88°, Δα = 4.22° и частотой f = 10.1 Гц, б) α₀ = 15.04°, Δα = 4.16° и частотой f = 10.05 Гц.

 

Положение точки отрыва потока можно извлечь из результатов CFD-расчета непосредственно по условию обнуления значения коэффициента трения (обнуления нормальной производной тангенциальной компоненты скорости) на поверхности профиля, что позволяет оценить правомерность применения формулы (2.1). На рис. 12 маркерами показана расчетная статическая зависимость для C_y.

 

Рис. 12. Сравнение расчетного значения x_s0 и C_y с решением Кирхгофа: 1 — CFD-расчет, 2 — аппроксимация по формуле Кирхгофа, 3 — аппроксимация по формуле Кирхгофа с x_s0 из расчета, 4 — x_s0 из расчета.

 

С использованием этих данных по формуле (2.1) для линейного случая вычислено положение точки отрыва x_s0 (нижний график, сплошная линия), а затем по этому положению снова вычислен C_y (верхний график, сплошная линия) по (2.1). Результаты для C_y весьма близки. На нижнем графике рис. 12 пунктирной линией показана полученная в расчете зависимость x_s0(α). Видно, что в нелинейном случае отрыв потока начинается позднее, при α = 11° по сравнению с α = 9° в линейном случае.

На верхнем графике рис. 12 пунктирной линией показано решение Кирхгофа для C_y, если в качестве x_s0 выбрана расчетная кривая. Она превышает расчетные значения C_y (маркеры) при α > 9°. Однако в целом согласование удовлетворительное, и можно заключить, что, во-первых, внутренняя переменная x_s действительно описывает положение точки отрыва потока, а во-вторых, что использование линейного приближения Кирхгофа для C_y(α, x_s) оправдано.

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

$\begin{array}{l}{C}_y\left(t\right)=C_{y_0}+C_{{}_y}^{\alpha }\left(\alpha \right)\alpha +C_{{}_y}^{\overline{\dot{\alpha }}}\left(\alpha \right)\overline{\dot{\alpha }},\\ m_z\left(t\right)=m_{z_0}+m_{{}_z}^{\alpha }\left(\alpha \right)\alpha +m_z^{\overline{\dot{\alpha }}}\left(\alpha \right)\overline{\dot{\alpha }}.\end{array}$ (2.4)

Представляет интерес, как соотносятся результаты нелинейного моделирования (2.3) и традиционного подхода (2.4). Для сравнения было задано гармоническое движение $\alpha ={\alpha }_0+\Delta \alpha \sin \left(\omega t\right)$ с амплитудой Δα = 3° около средних углов атаки в диапазоне α = 0–20° с шагом 2° с частотами f = 0.5, 1, 2 Гц. Результаты моделирования с использованием выражения (2.3) показаны на рис. 13. По эти данным с помощью метода линейной регрессии были идентифицированы коэффициенты модели (2.4).

 

Рис. 13. Моделирование C_y, m_z при колебаниях с малой амплитудой (а), аэродинамические производные m_z (в), моделирование колебаний с малой амплитудой при помощи аэродинамических производных (б).

 

Для коэффициента момента тангажа статическая $m_z^{\alpha }$ и нестационарная $m_z^{\overline{\dot{\alpha }}}$ показаны на рис. 13в. Линии для f = 1 и 2 Гц эквидистантно сдвинуты по вертикали на 0.1 и 0.2 для ясности рисунка. Полученные кривые демонстрируют типичные особенности, наблюдаемые в эксперименте в аэродинамической трубе с аналогичным движением. А именно: демпфирование тангажа остается примерно постоянным при безотрывном обтекании, затем по мере развития отрыва потока его величина нелинейным образом уменьшается и восстанавливается до прежних значений при полностью отрывном течении. Минимальное значение $m_z^{\overline{\dot{\alpha }}}$ зависит от частоты колебаний, что указывает на нелинейный характер зависимости m_z от параметров движения и, следовательно, неприменимость модели (2.4) в этой области.

Подставляя полученные коэффициенты обратно в формулы (2.4), можно получить зависимости коэффициентов силы и момента от времени C_y(t), m_z(t), какими они были бы в случае справедливости линейного приближения (см. рис. 13б). Можно видеть, как и следовало ожидать, различия в области развития отрывного обтекания, в особенности для большой частоты f = 2 Гц.

ЗАКЛЮЧЕНИЕ

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

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

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

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

ФИНАНСИРОВАНИЕ

Работа выполнена при финансовой поддержке Российского научного фонда (проект № 21-19-00659).

КОНФЛИКТ ИНТЕРЕСОВ

Авторы работы заявляют, что у них нет конфликта интересов.

作者简介

К. Abramova

Zhukovski Central Aerohydrodynamic Institute (TsAGI)

编辑信件的主要联系方式.
Email: kseniya.abramova@tsagi.ru
俄罗斯联邦, Moscow region

D. Alieva

Zhukovski Central Aerohydrodynamic Institute (TsAGI)

Email: diana.alieva@tsagi.ru
俄罗斯联邦, Moscow region

V. Sudakov

Zhukovski Central Aerohydrodynamic Institute (TsAGI)

Email: vit_soudakov@tsagi.ru
俄罗斯联邦, Moscow region

А. Khrabrov

Zhukovski Central Aerohydrodynamic Institute (TsAGI)

Email: khrabrov@tsagi.ru
俄罗斯联邦, Moscow region

参考

补充文件