CC BY
630
УДК 532.516
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ОТРЫВНЫХ ТЕЧЕНИЙ НА ОСНОВЕ НЕСТАЦИОНАРНЫХ УРАВНЕНИЙ НАВЬЕ-СТОКСА Д.А. Редчиц
Институт транспортных систем и технологий НАНУ,
Днепропетровск, Украина
Для численного моделирования нестационарных турбулентных отрывных несжимаемых течений применяются осредненные по Рейнольдсу уравнения Навье-Стокса. При моделировании турбулентности используются однопараметрические дифференциальные модели турбулентности. Решение системы исходных уравнений получено с помощью неявного конечно-объемного численного алгоритма, базирующегося на методе искусственной сжимаемости и многоблочных вычислительных технологиях. Тестирование разработанных методик проведено на задачах о циркуляционном течении в квадратной каверне, обтекании неподвижного и вращающегося цилиндров. Проведено сравнение результатов расчетов обтекания неподвижного и колеблющегося профилей. Представлены результаты расчета роторов Дарье и Савониуса с различным количеством и геометрическими характеристиками лопастей. Выполнен анализ поля течения вокруг роторов. Выделены основные стадии формирования вихревой структуры. Установлено влияние числа Рейнольдса, коэффициентов быстроходности и заполнения на энергетические характеристики роторов Дарье и Савониуса.
Ключевые слова: отрывное течение, число Рейнольдса, уравнения Навье-Стокса, несжимаемая жидкость, вихри, турбулентность, роторы Дарье и Савониуса.
Введение
Интенсивное развитие авиации, кораблестроения, транспорта привлекло большое внимание к исследованию отрывных течений при различных значениях режимных параметров и конфигурациях обтекаемых тел. Отрывное течение - часто встречающийся и, вместе с тем, наиболее сложный для исследования вид движения реальной жидкости [7]. Главная особенность отрывных течений заключается в том, что после появления отрыва потока течение становится нестационарным [3,25,26]. Несущие способности крыльев, аэродинамические характеристики летательных аппаратов, кораблей и подводных лодок, эффективность гидромашин и ветроустановок непосредственно зависят от развития отрыва потока. Зарождение и развитие отрыва потока, характеристики отрывных течений определяются большим числом параметров. Изучению этого сложного явления посвящено много экспериментальных [25,
26, 41] и теоретических работ [6, 7, 27].
Первые теоретические результаты были получены с помощью классической теории пограничного слоя, которая, однако, не позволяет учитывать
сильное взаимодействие пограничного слоя с внешним потоком. Наиболее важные теоретические и практические результаты по расчету течений с отрывом потока были получены с помощью асимптотической теории вязкой жидкости и полуэмпирических методов, базирующихся на априорном выборе модели течения на основе численного или физического экспериментов [6]. Среди различных подходов, применяемых к решению этой проблемы, важное место занимает математическое моделирование на основе уравнений Навье-Стокса. В настоящее время его роль возрастает с развитием ЭВМ, совершенствованием используемых моделей и численных методов, а также в связи с возможностью заменить расчетом дорогостоящий, а в ряде случаев практически невозможный физический эксперимент. Дополняя друг друга, расчет и эксперимент предоставляют новые возможности для изучения сложных взаимозависимых процессов.
Основная проблема получения нестационарных решений уравнений Навье-Стокса несжимаемой жидкости заключается в трудностях одновременного решения уравнений количества движения и уравнения неразрывности. На первом этапе развития численных алгоритмов решения уравнений Навье-Стокса для несжимаемых течений чаще использовались переменные завихренность-функция тока [22, 9, 27]. На основе данного подхода было решено большое количество прикладных задач, но расчеты пространственных задач с использованием функций тока весьма затруднительны.
Применение физических переменных позволяет решать двумерные и трехмерные задачи по единому алгоритму. Основные математические проблемы при решении уравнения Навье-Стокса связаны с различным типом дифференциальных уравнений для законов сохранения массы и количества движения. Различные способы преодоления указанных трудностей связаны с использованием для определения давления специального уравнения Пуассона [1, 37], уравнений для поправок [1, 14], различных штрафных функций [24], дополнением уравнения неразрывности нестационарным членом [37], регуляризацией матрицы коэффициентов при производных по времени [4, 8, 12, 23,
29, 30].
При решении уравнений несжимаемой жидкости в физических переменных применялся метод маркеров и ячеек (МАС) [10, 34, 35, 37], алгоритм SIMPLE [14, 33, 42], метод искусственной сжимаемости [37], в последнее время значительный прогресс в повышении эффективности численных алгоритмов достигнут при использовании локального предобуславливания [49, 55]. Обзор указанных методов можно найти в работах [15, 16].
В последние годы для моделирования течения несжимаемой жидкости все чаще применяются подходы, использующие эффекты сжимаемости [36, 50]. Практика применения уравнений сжимаемого газа показывает низкую работоспособность данной модели при малых числах Маха. Это связано с жесткостью исходных уравнений для низкоскоростных течений вследствие значительных различий в характерных временах конвективного переноса и распространения акустических возмущений. Данное препятствие может быть преодолено применением предобуславливания (preconditioning) [36, 50].
Метод искусственной сжимаемости представляет собой разумный компромисс между указанными выше подходами. С одной стороны, за счет добавления к уравнению неразрывности производной давления по времени исходная система уравнений приводится к единому типу. Это позволяет напрямую согласовать поля давления и скорости на одном временном слое.
Выбор модели турбулентности для расчета отрывных течений составляет отдельную проблему. Несмотря на существенный прогресс в моделировании турбулентности крупными (LES) и отсоединенными (DES) вихрями, прямым численным моделированием (DNS), при решении практических задач широко используются только модели на основе осредненных по Рейнольдсу уравнений Навье-Стокса (RANS). При замыкании RANS используются алгебраические, одно- и двухпараметрические модели турбулентной вязкости. Большую популярность получила однопараметрическая модель Спаларта-Аллмараса, демонстрирующая разумный компромисс между вычислительными затратами и точностью получаемых результатов.
В настоящей работе для исследования нестационарных турбулентных несжимаемых отрывных течений применяются осредненные по Рейнольдсу уравнения Навье-Стокса. Основное внимание в работе уделено исследованию развития отрыва несжимаемого потока при обтекании неподвижных, колеблющихся и вращающихся цилиндра и профиля, роторов Дарье и Савониуса.
1. Постановка задачи
Исходные уравнения. Для описания нестационарных турбулентных отрывных несжимаемых течений воспользуемся осредненными по Рейнольдсу уравнениями Навье-Стокса несжимаемой жидкости
Й=0’ (1)
Bui B (Ujui)
^ dxj
1 Bp д
Bui duj \
dt ' dxj p dxi ^ dxj e^ ^ dxj ^ dx¡) ’
где Xi, i = 1, 2 - декартовые координаты; t - время; ui - декартовы составляющие вектора средней скорости; p - давление; veff = v + vt - эффективный коэффициент кинематической вязкости; v и vt - молекулярный и турбулентный коэффициенты кинематической вязкости.
Моделирование турбулентности. Для замыкания осредненных по Рейнольдсу уравнений Навье-Стокса используются дифференциальные однопараметрические модели Spalart-Allmaras [53], SARC [52, 54] и SALSA [48], которые разработаны для задач внешней дозвуковой аэродинамики.
Однопараметрическая модель турбулентности Spalart-Allmaras (SA) [53] разработана в 1992 году и предназначена для описания равновесных течений типа пограничного слоя для задач внешнего обтекания при малых углах атаки с небольшими отрывными зонами. Генерация турбулентности определяется ротором поля скорости.
В 1997 году Spalart и Shur предложили модификацию исходной модели SARC [52, 54] введением дополнительных эвристических функций для учета кривизны линий тока и вращения твердой обтекаемой поверхности.
В 2003 году сотрудниками TU Berlin была предложена модификация SALSA [48], в которой генерация турбулентности связывалась не с ротором поля скорости, а с тензором скоростей деформаций. Кроме того, на основе опыта использования первоначальной модели были внесены изменения в константы. Различного рода модификации позволили более точно учитывать нестационарные эффекты в турбулентности.
Стандартная модель турбулентности Спаларта-Аллмараса (Spalart-Allmaras - SA) [53] предназначена для определения размерного кинематического коэффициента турбулентной вязкости
Vt = Vt • f
vi
(3)
где ^1 = х3/ (х3 + с30 - демпфирующая функция кинематических вязкостей х = . Здесь - рабочая переменная. Уравнение для определения
в модели Спаларта-Аллмараса имеет вид [53]
Dvt _ г- , 1 9
— _ CbiSvt +
(v + Vt)
дщ
Bxk
0,2 dvt dvt ~T~ Ju
a BxkBxk
Cbl
k2
1 + cb2
a
’¿V
_d_
2
Первое слагаемое в правой части (4) - источниковый член генерации турбулентности, где
ее ,Д,3Ж + ^Гр/,.2, у? = \/ЩЩ, (5)
Wij = 0.5 (Зщ /дх^ — дп^ /дх{) - тензор завихренности. Функция /у2 определяется с помощью соотношения
/у2 = 1 — ^
1 + хМ
Второе и третье слагаемые в правой части (4) отвечают за диссипацию турбулентности. Четвертое - за деструкцию турбулентности вблизи твердой стенки и содержит функцию
fw g
1 + с6
w3
g6 + с
w3
1/6
(7)
где g = г + cw2{r& — г), г = -з—-—. Здесь d - расстояние до ближайшей стенки.
S k2d2
Значения констант: к = 0.41 - константа Кармана, a = 2/3 - турбулентное число Прандтля, cb1 = 0.1355, cb2 = 0.622, cv2 = 5.0, cw2 = 0.3, cw3 = 2,
cv1 7.1, fv3 1-
Модель турбулентности Спаларта-Аллмараса, с учетом вращения и кривизны линий тока (Spalart-Allmaras for Rotation and Curvature - SARC) разработана Spalart и Shur [52]. В модели SARC [52], в отличие от стандартной SA [53], слагаемое генерации турбулентности умножается на эвристическую функцию fr1,
2r*
fr 1 = (1 + Cri) ——- [1 - cr3 tan“1 (cr2f)] - cr 1, (8)
(1 + r*)
где cr1 = 1, cr2 = 12. Константа cr3 определена Spalart и Shur [52] в диапазоне
0.6 ^ 1.0. В настоящей работе cr3 = 0.6. Функция r* определяется как r* = S/W. Здесь
S = у/2SijSij, W = y/2WijWij, Sl3 = 0.5 (диг/дх3 + ди3/дхг)
- тензор скоростей деформаций. Функция r согласно [52] имеет вид
r = 2WlkS]k
DS
г3 I (r Q. i p. C. 'IQ
D4 , (9)
где Qm - значение угловой скорости вращения относительно оси xm, єітп -компоненты тензора третьего ранга Леви-Чивита, D = \/0.5 (S2 + W2).
Модель Спаларта-Аллмараса, адаптированная к тензору скоростей деформаций (Strain-Adaptive Linear Spalart-Allmaras Model - SALSA) [48] является развитием оригинальной модели SA. Она основана на принципе вихревой вязкости для слабосжимаемых течений с пренебрежимо малыми флуктуациями плотности.
Уравнение для определения в модели Спаларта-Аллмараса, адаптированной к тензору скоростей деформаций, записывается в виде [48]
Dvt
Dt
— CbiSvt +
д
dxk
* 1 сЪ2 дщ дщ — fw C-bí 1 + 0>2
V 1 а) дхк_ a dxkdxk к2 а _d_
(10)
От стандартной модели турбулентности Спаларта-Аллмараса SALSA отличается модификацией членов генерации, диссипации и деструкции турбулентной вязкости. Стандартный коэффициент вы модифицируется в слагаемом генерации следующим образом
2
сы = сыл/Г, Г = тій (1.25, max(7, 0.75)), 7 = max(ai, а2), (11)
\ 0.65 0 65
1'01шЬО ’ "2=max(°> i_tanh(<l)) ■ (12)
Основное отличие SALSA от стандартной модели турбулентности SA заключается в использовании тензора скоростей деформаций вместо тензора завихренности, а именно
S — S *
— \ + fv 1
х)
дгіі dui
_ 1 дщ Здхк У (13)
Функции, которые входят в слагаемое деструкции, имеют следующий вид
dxj dxi
r — 1.6 tanh
0.7
Ф'
~S
Ф —
(14)
Значения всех остальных констант такие же, как и в стандартной модели БЛ.
Начальные и граничные условия. В качестве начальных условий задавались параметры невозмущенного потока во всей расчетной области. На
внешней границе применялись неотражающие граничные условия. На поверхности твердого тела ставилось условие прилипания. В моделях турбулентности SA, SARC, SALSA значение рабочей переменной на теле задавалось равным нулю щ = 0, на входной границе = 0.1, на выходной - ставилось условие Неймана.
2. Программно-методическое обеспечение для численного моделирования нестационарных течений на основе уравнений Навье-Стокса
Разработка программно-методического обеспечения для численного моделирования нестационарных несжимаемых турбулентных течений на основе уравнений Навье-Стокса включает следующие этапы: запись исходных уравнений в криволинейных неортогональных координатах, получение их дискретного аналога, построение расчетных сеток, решение системы алгебраических уравнений, визуализация результатов расчетов.
Система исходных уравнений (1,2), замкнутая одной из дифференциальных моделей турбулентности, интегрировалась численно с использованием метода контрольного объема. Разработанный численный алгоритм базируется на трехслойной неявной схеме Rogers-Kwak [45, 46], имеющей второй порядок точности интегрирования по времени, третий порядок противопоточ-ной аппроксимации конвективных слагаемых, и второй порядок центральноразностной аппроксимации диффузионных членов
dD dD d ^ ^ \ д
Imlk + Ж = (Е - Еи) ~ Щ (F ~ 1 = ~к'
где R - вектор невязок уравнений,
(15)
Р
u
v
É-1
ви
ixp + uU + 6tu _ 6yP + vU + 6v _
ev
ПхР + uV + ntu ПуР + vV + ntv
1т = diag [0, 1, 1] .
Вязкие члены в криволинейной системе координат имеют вид
Ev =
v + Щ J Re
0
(6X + ) uí + (6xVx + 6УЦу)
_ (й + 6У) vC + (6xVx + 6УЧу)
un
vn
У + Щ Л1е
(£хПх + £уПу) Щ + {п1 + п‘2) Щ , (17)
(£хПх + £уПу) ч + {п1 + ПЙ Щ _
где
,/ = М = аеЛ^ 4
д (х,у) [_ Пх Пу
- якобиан преобразования координат; £ = —хт£х — ут£у, п - —хтПх — УтПу,
- контравариантные компоненты вектора скорости; Яе - число Рейнольдса.
Для создания дискретного аналога исходных уравнений использовались регулярные сетки. В неодносвязных областях применялись многоблочные вычислительные технологии, в которых размерность отдельных пересекающихся сеток (блоков) не связана между собой. Такой подход позволил выработать единую методологию расчета течений вязкой жидкости при обтекании тел сложной геометрической формы.
В моделях турбулентности для аппроксимации конвективных слагаемых применялась схема ТУБ с ограничителем потоков ¡БКАБ третьего порядка. Производные в вязких членах аппроксимировались центрально-разностной схемой второго порядка.
Алгоритм решения уравнений базируется на трехслойной неявной схеме с подытерациями по псевдовремени т второго порядка точности по физическому времени £
£Х - 3уп 5 £у - 35 ПХ ------ '^у£ 5 пу - 3
метрические коэффициенты;
и - £ХП + £уV, V - ПхП + ПуV
1.5!п+1’т — 2!п + 0.5!п—1^
(18)
где верхний индекс п обозначает момент времени Ї = пД£. Для решения уравнений (4) и выполнения уравнения неразрывности на слое п+1 вводится псев-довременной слой т. Уравнения решаются итеративно, так чтобы ип+1,т+1 и ^п+і,т+і приближались к значению скорости мп+1, Цп+1 на новом временном слое, а дивергенция скорости стремилась к нулю. Полученная блочноматричная система алгебраических уравнений решалась методом итераций Гаусса-Зейделя.
Рассмотренная выше постановка задачи и численная методика реализованы в виде универсального комплекса программ, разработанного по модульному принципу, позволяющих включать новые модели турбулентности, численные методики, геометрии обтекаемых поверхностей, визуализировать структуру течений.
а)
б)
в)
г)
3. Верификация и тестирование программно-методического обеспечения
Современные требования к достоверности получаемых численных результатов и надежности программно-методического обеспечения требуют тщательного тестирования и верификации разработанного комплекса программ.
о 0 .
Рис. 1. Контуры завихренности при обтекании цилиндра (Ке=200) для моментов времени: а) Ь = 0; б) Ь = 2; в) Ь = 4; г) Ь = 6.
В разделе представлены результаты тестирования разработанной методики, алгоритмов и комплекса программ на задачах о развитии течения в квадратной каверне и обтекании неподвижного цилиндра несжимаемой жидкостью Кшак [16].
Получены структура стационарного и нестационарного ламинарного течения в следе за цилиндром, коэффициенты лобового сопротивления, подъемной силы и числа Струхаля. Показано, что при числах Рейнольдса Яе< 5 обтекание цилиндра происходит без отрыва потока. Увеличение числа Рейнольдса приводит к возникновению и развитию отрыва потока от поверхности цилиндра. При числах Рейнольдса Яе< 45 реализуется стационарный режим обтекания, который характеризуется наличием в ближнем следе двух симметричных вихрей. Обтекание цилиндра при Яе> 45 сопровождается образованием вихревой дорожки Кармана с числом Струхаля, зависящим от числа Рейнольдса (рис. 1). Результаты исследований хорошо согласуются с известными расчетными и экспериментальными данными.
Рис. 2. Профили скорости при обтекании кругового цилиндра (Ие= 100) в зависимости от линейной скорости поверхности: (а) а = 0; (б) а = 1.0; (в) а = 2.0.
4. Обсуждение результатов расчетов
На базе разработанного программно-методического обеспечения выполнено численное моделирование нестационарных турбулентных отрывных несжимаемых течений: обтекание вращающегося цилиндра [18], докритическое и закритическое обтекание неподвижного и колеблющегося профилей [19,20], роторов вертикально-осевых ветроагрегатов Дарье и Савониуса [21,20,43,44].
а)
б)
4.1. Ламинарное обтекание вращающегося цилиндра (эффект Магнуса)
Параметры задачи выбраны такими, при которых течение было ламинарным. Результаты расчетов показали, что вращение приводит к ускорению течения на одной стороне цилиндра и замедлению на другой в зависимости от величины а (отношение линейной скорости вращения поверхности цилиндра к скорости невозмущенного потока). Профили скорости показаны на рис. 3. Цилиндр вращается против хода часовой стрелки.
а)
б)
в)
Рис. 3. Контуры завихренности при обтекании кругового цилиндра (Б,е= 100) в зависимости от линейной скорости поверхности: (а) а = 0; (б) а = 1.0; (в) а = 2.0.
Контуры завихренности для различных значений а в момент времени, соответствующий минимуму подъемной силы при Ке= 100, показаны на рис. 3. С ростом скорости вращения вихри в следе, сошедшие с верхней стороны цилиндра, становятся более крупными, чем вихри, которые сходят с противоположной стороны (рис. 3б). При а > 2.0 наблюдается наличие двух стационарных вихрей, присоединенных к цилиндру (рис. 3в). Вращение цилиндра приводит к возникновению поперечной силы, величина и знак которой зависит от значений а (рис. 4).
С увеличением величины а давление на верхней части цилиндра становится большим, чем на нижней, и появляется ненулевая, осредненная по времени, поперечная сила (сила Магнуса). Фазовые диаграммы (зависимость С^ от Св) для Яе= 100 приведены на рис. 6. Замкнутость фазовых диаграмм
свидетельствует о периодической структуре течения. При значении а=1.9 фазовая диаграмма переходит в точку.
Таким образом, при постоянной скорости обтекания на цилиндр действуют аэродинамические силы, переменные по величине, направлению и частоте, зависящие от значений текущей угловой скорости.
4.2. Докритическое и закритическое обтекание профиля
NACA 4412 турбулентным потоком
Было выполнено сравнение результатов расчетов обтекания профиля NACA 4412 при числе Re= 1.64 х 106 с использованием моделей турбулентности SA, SARC и SALSA.
Установлено, что для докритического режима обтекания профиля (слабый отрыв - 12°) выбор модели турбулентности не оказывает существенного влияния на результаты расчетов.
На закритическом режиме обтекания профиля (массивный отрыв - 18°) наибольшие размеры отрывной зоны наблюдаются при использовании модели турбулентности SALSA (рис.7). Модель турбулентности SA генерирует более наполненный профиль турбулентной вязкости в пограничном слое по сравнению с моделями SARC и SALSA (рис. 8). Использование модели SARC приводит к незначительному улучшению получаемых результатов по сравнению с SA. Модель турбулентности SALSA лучше, чем SA и SARC передает нестационарную двумерную структуру течения с развитым отрывом потока.
Рис. 4. Изменение коэффициента подъем- Рис. 5. Фазовая диаграмма зависимости ной силы Съ от времени £ (Ие= 100). подъемной силы от лобового сопротивления
для различных значений а (Ие= 100).
в)
Рис. 6. Линии тока, построенные по мгновенному полю скоростей при угле атаки 18°: а) SA; б) SARC; в) SALSA.
в)
Рис. 7. Изополосы турбулентной вязкости при угле атаки 18°: а) SA; б) SARC; в) SALSA.
Распределение коэффициента давления по поверхности профиля для углов атаки 12° и 18° приведены на рис. 9. До угла атаки 12° наблюдается хорошее совпадение результатов по всем исследуемым моделям турбулентности с экспериментальными данными. При углах атаки, больших 12°, применение модели SALSA приводит к заметному улучшению получаемых результатов.
До угла атаки 12° результаты, полученные с использованием моделей турбулентности SA, SARC и SALSA, для коэффициентов подъемной силы и лобового сопротивления близки между собой и хорошо совпадают с экспериментальными данными (рис.10). Использование модели турбулентности SALSA приводит к существенному улучшению получаемых результатов в сравнении с моделями SA и SARC при закритическом режиме обтекания.
а) б)
Рис. 8. Распределение коэффициента давления по поверхности профиля NACA 4412 для различных углов атаки: (а) 12°; (б) 18°.
♦— Wadcock - — SA - -Д - SARC - -0 - SALSA ■
/ / // // //
к Г’
0 2 4 6
10 12 14 16 18 20
а)
б)
Рис. 9. Зависимость коэффициентов подъемной силы (а) и лобового сопротивления (б)
от угла атаки для профиля NACA 4412.
4.3. Обтекание осциллирующего профиля NACA 0015
Расчеты обтекания осциллирующего профиля NACA 0015 проведены при числе Рейнольдса Re= 1.95 х 106 для трех режимов обтекания: а) слабый отрыв потока, соответствующий среднему углу атаки ао = 4°; б) докритическое обтекание профиля (развитый отрыв), соответствующее ао = 11°; в) закри-тическое обтекание профиля (массивный отрыв), соответствующее ао = 15°. Мгновенный угол атаки крыла определялся по закону а (t) = ао + а\ sin (ut). Амплитуда колебаний составляла а\ = 4.2°, а безразмерная частота k = шс/2Уж = 0.1. Основное течение в случае докритического обтекания профиля крыла при наличии развитого отрыва потока (ао = 11°) - стационарное, отрывная зона не превышает половины длины профиля, наблюдаются отдельные колебания в следе и в части отрывной зоны. Зависимости коэффициентов подъемной силы и лобового сопротивления профиля от угла атаки
при гармонических колебаниях представлены на рис. 10. Движение по кривым направлено по часовой стрелке (рис. 10). Результаты, полученные с использованием модели турбулентности SALSA, удовлетворительно совпадают с экспериментальными данными.
а)
б)
Рис. 10. Зависимости коэффициентов подъемной силы (а) и лобового сопротивления (б) от углового положения профиля крыла (ао = 11 ).
На этих рисунках нанесены в соответствии со следующими обозначениями:
---модель SA — - модель SALSA
(настоящая работа); (настоящая работа);
о - эксперимент о - эксперимент
(статический профиль) [51]; (колеблющийся профиль) [40];
△ - модель Baldwin-Barth ж - модель Wilcox (расчет [32]); к — ш (расчет [32]);
□ - модель SA (расчет [32]); • - модель SALSA (расчет [32]).
Структура течения при закритическом обтекании профиля и наличии массивного отрыва потока (а0 = 15°) характеризуется ярко выраженными нестационарными явлениями (рис. 12). Отрыв потока зарождается на подветренной стороне вблизи носика профиля, распадаясь затем на систему вихрей с различными скоростями движения.
Разработанная методика воспроизводит структуру нестационарного отрывного обтекания осциллирующего профиля NACA 0015. Показано преимущество модели турбулентности SALSA по сравнению с моделями SA, SARC при расчете течений с развитым двумерным нестационарным отрывом потока.
Рис. 11. Контуры завихренности, полученные по модели турбулентности SALSA: а) 19° б) 19° в) 17° г) 14° | (|- профиль движется вверх, | - профиль движется вниз).
4.4. Роторы вертикально-осевых ветроустановок
Разработка и усовершенствование альтернативных источников энергии является актуальной проблемой энергетики. К одному из перспективных направлений решения данной проблемы относится ветроэнергетика. Большое распространение в мире получили двух- и трехлопастные горизонтальноосевые (ГО) ветроэнергетические установки (ВЭУ) пропеллерного типа. Это связано с высоким коэффициентом использования энергии ветра. Близкими значениями коэффициента мощности из вертикально-осевых (ВО) ВЭУ обладают только роторы Дарье.
Повышение мощности ВЭУ и увеличение коэффициента использования энергии ветра делает задачу выбора рациональной аэродинамической формы ротора весьма актуальной. Ведущую роль в работе ВЭУ играют нестационарные аэродинамические процессы, поэтому основным направлением исследований должна быть разработка новых универсальных методов расчета нестационарных процессов при обтекании потоком роторов ветроагрегатов.
Известные методики определения аэродинамических и энергетических характеристик роторов ВЭУ основаны на экспериментальных данных, импульсной и вихревой теориях, численном решении уравнения потенциала [5, 11, 41]. Они используют определенные допущения при постановке задачи (квазистационарность потока, отсутствие учета вязко-невязкого взаимодействия и т.д.). Главными трудностями в расчете нестационарных процессов при обтекании роторов ВО ВЭУ являются эффекты динамического срыва потока. До настоящего времени ни одна из известных упрощенных моделей не давала возможности адекватно рассчитать аэродинамические характеристики роторов в случае динамического срыва потока.
Уравнения Навье-Стокса являются одной из наиболее полных математических моделей механики жидкости и газа. Применение их совместно с дифференциальной моделью турбулентности, уравнением динамики ротора позволяет исследовать особенности нестационарного обтекания, структуру поля скоростей, динамический срыв потока, процессы формирования и распада вихрей вокруг самого ротора, а также в следе за ветроагрегатом. Одним из наиболее перспективных направлений расчета аэродинамических и энергетических характеристик ротора ВЭУ является совместное численное решение уравнений динамики вязкой несжимаемой жидкости и твердого тела.
Рассматрим ортогональные роторы Дарье и Савониуса, лопасти которых имеют длину, многократно превышающую хорду (рис. 12). В таком случае можно пренебречь концевыми эффектами на лопастях и воспользоваться гипотезой о плоскопараллельной структуре течений. Таким образом, задача обтекания ВО ВЭУ допускает двумерную постановку в плоскости, перпендикулярной оси вращения ротора. Роторы Дарье и Савониуса полагаем абсолютно твердыми. Поскольку для максимальных скоростей ветра и значений коэффициента быстроходности локальные числа Маха низкие (М < 0.3), поле течения принято несжимаемым.
Уравнение вращения ротора вертикально-осевой установки относительно неподвижной оси имеет следующий вид
^г~(й = ^ ~ ~ ^
где 1г - момент инерции ротора; и - угловая скорость вращения; Q - крутящий момент, обусловленный действием потока на лопасти ВЭУ; Qld - момент полезной нагрузки, приложенный к валу электрогенератора; Qfr - результирующий момент трения в электромеханической системе ветроагрегата.
Ротор Дарье. Ниже представлены результаты численного моделирования обтекания вращающейся одиночной лопасти, а также ротора Дарье с двумя и тремя лопастями.
Численное моделирование обтекания вращающейся лопасти проведено при различных коэффициентах быстроходности Л. По коэффициенту тангенциальной силы при Лі = 2.5 наблюдается широкий разброс экспериментальных и расчетных данных (рис. 13). Результаты работы лучше согласуются с известными экспериментальными данными, чем расчеты других авторов, особенно в наветренной части траектории лопасти.
а) б) Рис. 12. Расчетные схемы для роторов Дарье (а) и Савониуса (б).
Основной крутящий момент создается на наветренном участке траектории лопасти (рис. 14). Анализ результатов расчетов показал, что поток, проходя через наветренный участок траектории лопасти, теряет часть своей кинетической энергии. Именно поэтому коэффициент крутящего момента лопасти больше на этом участке, чем на подветренном. На подветренном участке траектории он минимальный (Ах = 2.5, А2 = 5.0) или, вообще, отрицательный (Аз = 7.5).
: -Л-- I
■А Ь * д :, * ЛД':/Л • л К' ;
<*' X. ^ г. 1 рЧ\; %
120
180
9
240
300
360
Рис. 13. Изменение коэффициента танген- Рис. 14: Изменение неосредненных коэф-
циальной силы лопасти Ст от углового фициентов крутящего момента Сд от
положения ротора 0 углового положения ротора 0 для различных
— - расчет [39]; о, +— эксперимент [39]; коэффициентов быстроходности Л.
△ - расчет [41];-----настоящая работа.
Для иллюстрации особенностей обтекания ротора Дарье были выбраны геометрические параметры и коэффициент быстроходности, соответствующие экспериментальной работе С. БгоеЫег [31] (рис. 15). На рис. 15б кроме стандартной визуализации вихрей добавлены сплошные и прерывистые линии, а также отдельные точки для того, чтобы стиль интерпретации расчетных данных соответствовал стилю визуализации экспериментальных данных работы С. БгоеЫег [31]. Приведена реконструкция структуры течения при работе двух- и трехлопастного роторов Дарье для коэффициента быстроходности А = 2.14 на основе натурного (а) и вычислительного (б) экспериментов (рис. 16). Для наглядности оставлены вихри максимальной интенсивности. Выделены стадии зарождения, развития и срыва вихрей при различных положениях лопасти на траектории.
В целом картина течения вблизи ротора Дарье характеризуется существенными нестационарными явлениями. К ним относятся: динамический срыв потока, образование сложной системы вихрей, повышение уровня турбулентности в затененной области, взаимодействие вихрей различных размеров, скоростей движения и интенсивности с подвижными поверхностями роторов. Полученная картина течения хорошо согласуется с имеющимися экспериментальными данными [31].
Установлено влияние чисел Рейнольдса, коэффициентов быстроходности и заполнения на энергетические характеристики ротора Дарье (рис. 17, 18). Показано, что рост числа Рейнольдса приводит к увеличению значений коэффициента мощности (рис. 17, 18). При уменьшении коэффициента заполнения ротора Дарье коэффициент мощности становится менее чувствительным к изменению коэффициента быстроходности.
Ротор Савониуса. При численном моделировании обтекании двух- и трехлопастных роторов Савониуса выполнены три типа вычислительных экспериментов. Первый тип - вычислительные эксперименты для неподвижного ротора Савониуса, который фиксировался при различных углах относительно набегающего потока с шагом АО = 10. Для большинства угловых положений ротора Савониуса осредненный по времени коэффициент крутящего момента положительный.
а)
а)
б) б)
Рис. 15. Визуализация течения при работе Рис. 16. Реконструкция структуры течения двухлопастного ротора Дарье для коэффи- при работе двухлопастного ротора Дарье для циента быстроходности Л = 2.14 на основе коэффициента быстроходности Л = 2.14 на натурного (а) и вычислительного (б) экспе- основе натурного (а) и вычислительного (б) риментов. экспериментов.
Второй тип - вычислительные эксперименты при фиксированном коэффициенте быстроходности ротора. Коэффициенты крутящего момента и мощности двух- и трехлопастного ротора Савониуса определялись осреднением за один полный оборот. Вращение ротора при Л =1.4 характеризуется квазиста-ционарным режимом обтекания. Выделены основные стадии формирования вихревой структуры при вращении ротора (рис. 19, 20). Периодичность в структуре течения вокруг ротора наблюдается через 180 и 120 для двухлопастного и трехлопастного соответственно. Визуализация обтекания выполнена с помощью контуров завихренности. Определены зависимости коэффициентов крутящего момента и мощности от коэффициента быстроходности.
У двухлопастного ротора значения энергетических характеристик выше, чем у трехлопастного (рис. 21). Полученные результаты удовлетворительно согласуются с известными экспериментальными данными [28].
0,4
- 0
Ср
-0,2
-0,4
« Яе = 104 П Во-1П!
д Ке = 106 : \
0,6
0,4
0,2
Ср
-0,2
-0,4
-0,6
о СГ=0.33 □ и = 0.67 д <7 = 1.00
Рис.17. Зависимость осредненного коэффициента мощности ротора Дарье от коэффициента быстроходности для различных чисел Рейнольдса (а = 0.67).
Рис. 18. Зависимость осредненного коэффициента мощности ротора Дарье от коэффициента быстроходности для различных коэффициентов заполнения (Кв= 105).
Третий тип вычислительных экспериментов - решение связанной задачи динамики и аэродинамики трехлопастного ротора Савониуса. Проанализирована картина течения вокруг ротора, приведены зависимости коэффициентов лобового сопротивления, подъемной силы и крутящего момента, а также угловой скорости вращения от времени (рис. 22).
Расчет проводился в три этапа. Целью первого этапа (Ь = 0 ^ 7) было получение периодического течения, по структуре похожего на дорожку Кармана. На втором (Ь = 7 ^ 13) и третьем (Ь = 13 ^ 23) этапах совместно с аэродинамической задачей решалось уравнение вращения ротора Савониуса. С момента времени Ь = 7 ротор освобождается и вращается под действием набегающего потока ветра.
Вращение ротора приводит к увеличению интенсивности вихрей. Частота схода вихрей определяется скоростью набегающего потока, характерными размерами и частотой вращения самого ротора.
На третьем этапе, в момент времени Ь = 13, к ротору Савониуса прикладывается момент нагрузки. Происходит стабилизация угловой скорости вращения ротора (относительно среднего значения и = 2.8), а также возникают близкие к периодическим, колебания коэффициентов лобового сопротивления, подъемной силы и крутящего момента (рис. 22).
Рис. 19: Контуры завихренности возле Рис. 20. Контуры завихренности возле
подвижного (Л = 1.4) двухлопастного ро- подвижного (Л = 1.4) трехлопастного ротора Савониуса: а) в = 0°; б) в = 45°; тора Савониуса: а) в = 0°; б) в = 30°;
в) в = 90°; г) в = 135°; д) в = 180°. в) в = 60°; г) в = 90°; д) в = 120°.
0,3
0,2 -
Ср
0,1
d □ 8 й' ft Ш 1 4] □
д Л □ :
0 0,2 0,4 0,6 д 0,8 1 1,2 1-4 а) 0 °-2 °-4 °-6 Я 0,8 1 1,2 1,4 б)
Рис. 21. Зависимость осредиениого за один оборот коэффициента мощности Ср от коэффициента быстроходности Л двух-(а) и трехлопастного (б) ротора Савониуса □ - эксперимент B. Blackwell (Re= 4.32 х 105, Re= 8.64 х 105); А, ■ - настоящая работа (Re= 4.32 х 105,
Re= 8.64 х 105).
Рис.22: Изменение неосредненных коэффициентов лобового сопротивления С0, подъемной си-
лы Сь, крутящего момента С<^ и угловой скорости вращения ш трехлопастного ротора Савониуса
от времени £ и углового положения ротора в.
5. Выводы
Для численного моделирования нестационарных турбулентных отрывных несжимаемых течений применяются осредненные по Рейнольдсу уравнения Навье-Стокса. Замыкание системы уравнений осуществлено с помощью модели турбулентности Спаларта-Аллмараса и ее модификаций. Реализация используемого подхода выполнена с помощью разработанного программнометодического обеспечения численного решения уравнений Навье-Стокса несжимаемой жидкости в произвольных неортогональных координатах на по-
движных сетках. Выполнено тестирование программно-методического обеспечения на задачах о циркуляционном течении в квадратной каверне, обтекании неподвижного и вращающегося цилиндров. Проведено сравнение результатов расчетов обтекания неподвижного и колеблющегося профилей. Представлены результаты расчета роторов Дарье и Савониуса с различным количеством и геометрическими характеристиками лопастей. Выполнен анализ поля течения вокруг роторов. Выделены основные стадии формирования вихревой структуры. Вязкие и динамические эффекты играют основную роль в работе ротора Дарье, максимальный крутящий момент создается на наветренном участке траектории лопасти. Установлено влияние числа Рейнольдса, коэффициентов быстроходности и заполнения на энергетические характеристики ротора Дарье.
Литература
1. Белов И.А. Задачи и методы расчета отрывных течений несжимаемой жидкости / И.А.Белов, С.А.Исаев, В.А.Коробков. - Л.: Судостроение, 1989.- 256с.
2. Белов И.А. Теплоотдача и сопротивление пакетов труб / И.А.Белов, Н.А.Кудрявцев. - Л.: Энергоатомиздат, 1987. - 223с.
3. Белоцерковский С.М. Математическое моделирование плоскопараллельного обтекания тел /С.М.Белоцерковский, В.Н.Котовский, М.И.Ништ. -М.: Наука, 1988. - 232с.
4. Блинова Л.А., Шур М.Л. Метод "масштабирования сжимаемости" для расчета нестационарных течений вязкого газа в широком диапазоне изменения характерных чисел Маха / Конструирование алгоритмов и решение задач математической физики / Л.А.Блинова. - М.: ИПМ, 1991.
- С.34-39.
5. Волков Н.И. Аэродинамика ортогональных ветродвигателей (некоторые математические модели и численная реализация). Учебное пособие /
Н.И.Волков. - Сумы: ВВП "Мрия-1" ЛТД, 1996. - 198с.
6. Гогиш А.В., Нейланд В.Я., Степанов Г.Ю. Теория двумерных отрывных
течений // Итоги науки и техники. Гидромеханика. - М.: Наука, 1975. -
8. - С.5-73.
7. Гогиш Л.В. Турбулентные отрывные течения / Л.В.Гогиш, Г.Ю.Степанов.
- М.: Наука, 1979. - 368с.
8. Гончаров В.А., Кривцов В.М., Чарахчьян А.А. Численная схема моделирования дозвуковых течений вязкого сжимаемого // ЖВМ и МФ. -1988. - 28;12. - С.1858-1866.
9. Госмен А.М. Численные методы исследования течений вязкой жидкости / А.М.Госмен, Пан В.М., Ранчел А.К. и др. - М.: Мир, 1972. - 323с.
10. Госмен А.М. Дальнейшее развитие метода маркеров и ячеек для течений несжимаемой жидкости / Численные методы в механике жидкостей. -М.: Мир, 1973. - С.165-173.
11. Кривцов В.С. Неисчерпаемая энергия / В.С.Кривцов, А.М.Олейников, А.И.Яковлев. - Харьков: Нац. аэрокосм. ун-т "ХАИ", 2003. - 919с.
12. Лапин Ю.В. Внутренние течения газовых смесей / Ю.В. Лапин, М.Х.Стрелец. - М.: Наука, 1989. - 368с.
13. Лойцянский Л.Г. Механика жидкости и газа /Л.Г.Лойцянский. - М.: Наука. 1987.- 840с.
14. Патанкар С.В. Численные методы решения задач теплообмена и динамики жидкости / С.В.Патанкар. - М.: Энергоатомиздат, 1984. - 152с.
15. Приходько А.А. Компьютерные технологии в аэрогидродинамике и тепломассообмене / А.А.Приходько. - Киев: Наукова думка, 2003. - 380с.
16. Приходько A.A., Редчиц Д.А. Численное моделирование нестационарного течения в следе за цилиндром на основе уравнений Навье-Стокса // Прикладная гидромеханика. - 2005. - 7;1. - С.56-71.
17. Приходько А.А., Редчиц Д.А. Математическое моделирование динамики и аэродинамики вертикально-осевых ветроагрегатов // Вестник Харьковского национального университета. - 2005. - 703;5. - С. 178-197.
18. Приходько А.А., Редчиц Д.А. Численное моделирование эффекта Магнуса на основе уравнений Навье-Стокса // Вісник Дніпропетровського університету. Механіка.- 2005. -1;7. - С.40-60.
19. Приходько A.A., Редчиц Д.А. Численное моделирование дозвукового обтекания осциллирующего профиля на основе уравнений Навье-Стокса // Техническая механика.—2006. - 1.- С.104—114.
20. Приходько A.A., Редчиц Д.А. Компьютерное моделирование аэродинамики подвижных роторов ветроагрегатов Дарье и Савониуса // Аэрогидродинамика: проблемы и перспективы. — 2006. — 2. — С. 120—142.
21. Редчиц Д.А., Приходько А.А. Численное решение связанной задачи динамики и аэродинамики ротора ветроагрегатов // Космическая наука и технология. — 2005. — 11;1. — С.27—35.
22. Роуч П. Вычислительная гидродинамика / П.Роуч. — М.: Мир, 1980.— 616с.
23. Стрелец М.Х., Шур М.Л. Метод масштабирования сжимаемости для расчета стационарных течений вязкого газа при произвольных числах Маха // ЖВМ и МФ. — 1988. — 28; 2. — С.254—266.
24. Темам Р. Уравнения Навье-Стокса. Теория и численный анализ / Р.Те-мам. — М.: Мир, 1981. — 408с.
25. Чжен П. Управление отрывом потока. Экономичность, эффективность, безопасность / П.Чжен. — М.: Мир, 1979. — 552 с.
26. Чжен П. Отрывные течения / Чжен П. — М.: Мир, Т.1, 1972 — 300с.; Т.2, 1973. — 280 с.; Т.3,1973. — 334с.
27. Шлихтинг Г. Теория пограничного слоя / Г.Шлихтинг. — М.: Наука, 1974.
— 711c.
28. Blackwell B.F., Sheldahl R.E., Feltz L.V. Wind tunnel performance data for two- and three-bucket Savonius Rotors // Sandia National Laboratories Albu-querque. SAND76-0131. — 1976. — P.105.
29. Briley W.R., McDonald H. Solution of the multidimensional Navier-Stokes equations by a generalized implicit method // Journal of Computation Physics. —1977. — 24;4. — P.372-397.
30. Briley W.R., McDonald H., Shamreth S.J. A low Mach number Euler formulation and application to time-iterative LBI schemes // AIAA Journal. — 1983. — 21;10. — P.1464—1469.
31. Brochier G., Fraunie P., Beguier C., Paraschivoiu I. Water channel experiments of dynamic stall on Darrieus wind turbine blades // Journal Propulsion.
— 1986. — 2;5. — P.445-449.
32. Bunge U., Martin A., Schmidt S., Schatz M., Thiele F. DES and its Applications at Technical University of Berlin // Proc. International Conference on DES. — Workshorp. — St. Petersburg, 2003.
33. Caretto L.S., Gosman A.D., Patankar S.V. and Spalding D.B. Two calculation procedures for steady three-dimensional flows with recirculation // Proceedings of the 3rd Int. Conference on Numerical Methods in Fluid Dynamics, Paris, France, 1972. — P. 60.
34. Chan R.K.-С., Street R. Lt., Strelkoff T. Computer studies of finite-amplitude water waves // Tech. Rep, 104, Dep. of Civil Eng„ Stanford Univ., Stanford, California.—1969.— Р.126.
35. Chan R.K.-С., Street R.L. A computer study of finite-amplitude water waves // Journal of Computation Physics. — 1970. — 6. — P.68-94.
36. Choi Y.-H., Merkle, C.L. The Application of Preconditioning to Viscous Flows // Journal of Computational Physics. — 1993. — 105. — P.207-223.
37. Chorin A.J. A numerical method for solving incompressible viscous flow problems // Journal of Computation Physics. — 1967. — 2. — P.12-26.
38. Harlow F.H., Welch J.E. Numerical calculation of time-dependent viscous incopressible flow with free surface // Physics of Fluids. — 1965. — 8;12. — P.2182—2189.
39. Oler J.W., Strickland J.H., Im B.J., Graham G.H. Dynamic stall regulation of the Darrieus turbine // SAND83-7029. Texas technical university. — 1983.
— P.154.
40. Piziali R.A. An experimental investigation of 2D and 3D oscillating wing aerodynamics for a range of angle of attack including stall // NASA TM 4632. —1993.
41. Paraschivoiu I. Wind turbine design with emphasis on Darrieus concept /
I.Paraschivoiu. — Canada: Polytechnic international press, 2002. — 438p.
42. Patankar S.V., Spalding D.B. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows // Journal Heat and Mass Transfer. - 1972. - 15. - P.1787-1806.
43. Prykhodko O., Redchyts D. Mathematical modelling of vertically-axis wind turbine work on the basis of Navier-Stokes equations // Proc. Flow and Transport Processes in Complex Obstructed Geometries: From Cities and Vegetative Canopies to Industrial Problems. - Institute of Hydromechanics of NASU Kyiv (Ukraine). - 2004. - P.161.
44. Prykhodko O., Redchyts D. Numerical modelling of dynamics and aerodynamics processes of Darrieus and Savonius rotors // Proc. 77th Annual Meeting of the Gesellschaft fur Angewandte Mathematik und Mechanic. - Technische Universitat Berlin. - 2006. - P.346.
45. Rogers S.E., Kwak D. An upwind differencing scheme for the incompressible Navier-Stokes equations // Journal Numerical Mathematics. - 1991.- 8. -P.43-64.
46. Rogers S.E., Kwak D. An upwind differencing scheme for the time-accurate incompressible Navier-Stokes equations // AIAA Journal. - 1990. - 28;2. -P.253-262.
47. Roshko A. On the development of turbulent wakes from vortex streets // NACA Report. - 1954. - 1191. - P.32-65.
48. Rung T., Bunge U., Schatz M., Thiele F. Restatement of the Spalart-Allmaras eddy-viscosity model in strain-adaptive formulation // AIAA Journal. -2003. - 4;7. - P.1396-1399.
49. Turkel E. Review of preconditioning methods for fluid dynamics // Applied Numerical Mathematics. - 1993. - 24. - P.257-284.
50. Turkel E., Vatsa V. N., Radespiel R. Preconditioning Methods for Low-Speed Flows // ICASE. - 1996. - 57. - 19pp.
51. Sheldahl R.E., Climas P.C. Aerodynamic characteristics of seven symmetrical airfoil sections through 180-degree angle of attack for use in aerodynamic analysis of vertical axes wind turbines // Sandia National Laboratories Albuquerque. SAND80-2114. -1995. - P.118.
52. Shur M.L, Strelets M. K., Travin A.K., Spalart P.R. Turbulence modeling in rotating and curved channels: Assessing the Spalart-Shur correction // AIAA Journal. - 2000. - 38;5. - P.784-792.
53. Spalart P.R., Allmaras S.R. A one-equation turbulence model for aerodynamic flow // AIAA Paper. - 1992.- 12;1. - P.439-478.
54. Spalart P.R., Shur M. On the sensitization of turbulence models to rotation and curvature // Aerospace science and technology Journal. - 1997. - 1;5. -P.297-366.
55. Van Leer B., Lee W.-T., Roe P. Characteristic timestepping or local preconditioning of the Euler equations // AIAA Paper. - 1991. - 19;6. - P.1552-1587.
MATHEMATICAL MODELLING OF SEPARATED FLOWS ON THE BASIS OF UNSTEADY NAVIER-STOKES EQUATIONS D.A. Redchits
Institute of Transport Systems and Technologies NANU,
Dnepropetrovsk, Ukraine
The Reynolds averaged Navier-Stokes equation is applied to the numerical simulation of unsteady turbulent separated incompressible flows. The one-parameter differential models of turbulence models are used. The solution of the initial equation system is obtained using the implicit finite-volume numerical algorithm which is based on the artificial compressibility method. Besides, overlapped structured grids are applied in multidomain areas. Verification of developed CFD code is fulfilled on problems concerning the circulating flow in the square cavity, the flow near fixed and rotated cylinders. Comparison of of calculations concerning the flow near fixed and oscillating airfoil is carried out. The results of calculations of Darrieus and Savonius rotors with different number of blades and blades with differnt geometrical characteristics are proposed. The analysis of flow near rotors is done. Main stages of vortex development are allocated. The Reynolds number influence on power characteristics of Darrieus and Savonius rotors is found as well as tip-speed-ratio and solidity.
Key words: separated flow, Reynolds number, Navier-Stokes equations, imcompressible liquid, vortex, turbulence, Darrieur and Savonius rotors.
