CC BY
54
Методика решения задач аэроупругости для лопасти ветроустановки с использованием СПО
1 2 П.С. Лукашин <ski1l@mail.ru> 1 2 В.Г. Мельникова <vg-melnikova@yandex.ru> 2 С.В. Стрижак <strijhak@yandex.ru> 1 Г.А. Щеглов <shcheglov_ga@bmstu.ru> 1МГТУ им. Н.Э. Баумана, 105005, г. Москва, ул. 2-я Бауманская, д.5. стр.1 2 Институт системного программирования им. В.П. Иванникова РАН, 109004, Россия, г. Москва, ул. А. Солженицына, д. 25
Аннотация. В связи с развитием ветроэнергетики и строительством новых ветропарков в РФ возникает потребность в решении ряда прикладных задач и разработке эффективных методик расчета элементов ветроустановок. Одно из направлений в задачах механики сплошной среды связано с задачами аэроупругости. В данной статье показана возможность решения связанной задачи аэроупругости с использованием программного комплекса на базе свободного программного обеспечения OpenFOAM и Code_Aster. На примере лопасти ветроустановки длиной 61.5 метра рассмотрены методики решения задач статической и динамической аэроупругости, в которых расчет обтекания лопасти дозвуковым набегающим потоком воздуха производится в пакете OpenFOAM (решатели simpleFOAM и pimpleFOAM), а расчет напряженно-деформированного состояния лопасти производится в пакете Code_Aster. В статье приводятся блок-схемы для трех различных подходов решения задачи аэроупругости, примеры скриптов и командных файлов для передачи данных между пакетами в процессе расчета. Контрольно-объемная сетка, состоящая их гексаэдральных элементов, для расчета обтекания лопасти построена в пакете OpenFOAM, конечно-элементная сетка, состоящая из треугольных оболочечных элементов первого порядка, для расчета напряженно-деформированного состояния построена в пакете Salome-Meca. Результаты расчета представлены в форме полей для давления и скорости; зависимостей для невязкок давления, скорости, турбулентной вязкости; проекций аэродинамической силы от времени; эпюр перемещения и напряжения; значений давления и перемещения для двух точек на поверхности лопасти от времени. Расчеты выполнены с использованием ресурсов web-лаборатории UniHUB ИСП РАН.
Ключевые слова: аэроупругость; ветроустановка; лопасть; Code_Aster; OpenFOAM; модель турбулентности; решатель; статический расчет; динамический расчет; частота, перемещения, напряжение, расчеты, web-лаборатория.
DOI: 10.15514/ISPRAS-2017-29(6)-16
Для цитирования: Лукашин П.С., Мельникова В.Г., Стрижак С.В., Щеглов Г.А. Методика решения задач аэроупругости для лопасти ветроустановки с использованием СПО. Труды ИСП РАН, том 29, вып. 6, 2017 г., стр. 253-270. DOI: 10.15514/ISPRAS-2017-29(6)-16
1. Введение
В связи с развитием ветроэнергетики в РФ, проектированием новых ветроэлектрических установок (ВЭУ), ветропарков и их эксплуатацией в различных климатических условиях на обширной территории РФ возникает потребность в решении ряда прикладных задач. К одной из них можно отнести оценку динамических и прочностных характеристик лопастей ВЭУ с учетом ветровой нагрузки. Наиболее приспособленными для решения таких задач являются универсальные коммерческие программные пакеты. Но статистика показывает, что различные крупные и средние предприятия, постепенно стали осознавать преимущества, которые предоставляет свободное программное обеспечение (СПО): снижение затрат на программные комплексы, открытость исходного кода, а также защита от возможных санкций иностранных правообладателей. Поэтому поиск, модификация, усовершенствование существующих и создание новых программных средств с открытым исходным кодом является актуальной задачей на сегодняшний день.
В настоящее время СПО уже успешно используется для расчета конструкций в области энергетики. В частности, можно отметить пакет Code_Aster, используемый крупнейшей энергетической компанией Франции EDF. Важным преимуществом свободного программного обеспечения является возможность создания на базе нескольких пакетов нового программного комплекса для решения сложных задач мультифизики. К таким задачам традиционно относятся задачи аэроупругости.
Основной целью работы является отладка методики решения сопряженных задач взаимодействия деформируемой конструкции и потока среды на основе двух открытых пакетов: OpenFOAM и Code_Aster.
Возможны различные постановки задачи аэроупругости, основанные на упрощениях и допущениях [1,2]. Исходя из анализа литературы [2-5], можно выделить два варианта связанной задачи FSI:
• задачи статической аэроупругости: совместное рассмотрение уравнений теории упругости и уравнений для стационарного течения жидкости и газа;
• задачи динамической аэроупругости: совместное рассмотрение уравнений теории упругости и уравнений для нестационарного течения жидкости и газа.
В [6] выделяется несколько подходов к решению задачи FSI:
• one-way coupling: раздельное последовательное интегрирование уравнение двух подсистем на каждом шаге интегрирования, без учета перемещения тела в потоке;
• two-way coupling: согласованное интегрирование уравнений двух подсистем на каждом шаге интегрирования, с учетом перемещения тела в потоке.
Основываясь на характере действующей на конструкцию лопасти ВЭУ нагрузки, величине влияния деформации тела на поток и представленной выше классификации было выбрано три варианта методики расчета:
• квазистатическая постановка - расчет напряженно-деформированного состояния (НДС) лопасти при статическом нагружении полем давлений, полученным в результате стационарного аэродинамического расчета недеформированной конструкции.
• динамическая постановка - расчет вынужденных колебаний лопасти при нагружении полем давлений, полученным в результате нестационарного аэродинамического расчета.
• сопряженная квазистатическая постановка - расчет НДС лопасти при статическом нагружении полем давлений, полученным в результате стационарного аэродинамического расчета с учетом влияния деформации лопасти на параметры течения.
2. Описание методик
2.1. Квазистатическая постановка
Простейшую последовательную методику решения задачи аэроупругости, блок-схема которой показана на рис. 1, допустимо использовать в случае, когда режим обтекания лопасти можно считать безотрывным и стационарным, а влиянием упругих перемещений конструкции на характер течения и распределением аэродинамических нагрузок по лопасти можно пренебречь. Также упрощающим предположением является допущение о компенсации сил инерции диссипативными силами что приводит задачу определения НДС к статической. Методика состоит из следующих этапов:
• проведение стационарного аэродинамического расчета методом контрольного объема в пакете OpenFOAM с использованием решателя simpleFOAM. В результате определяются поля скоростей и давлений;
• перенос результатов аэродинамического расчета из пакета OpenFOAM в пакет Code_Aster;
• определение НДС конструкции методом конечных элементов в Code_Aster.
Здесь ключевым для совместного использования пакетов является второй этап, в ходе которого осуществляется проецирование поля давления, полученного в аэродинамическом расчете с контрольно-объемной (КО) сетки на конечно-элементную (КЭ), для дальнейшего использования в качестве нагрузки при расчете НДС.
Рис. 1. Алгоритм расчета при квазистатической постановке задачи Fig. 1. Simulation algorithm at quasistatic problem definition
Преобразование данных из формата *.foam, используемого в OpenFOAM, в формат *.med, используемый в качестве входного формата для Code_Aster осуществлялось посредством модуля визуализации Paraview, встроенного в препостпроцессор Salome-Meca. Для проецирования значений давлений с одной сетки на другую применялся оператор Code_Aster PROJ_CHAMP, в котором использовался метод коллокаций. В этом методе определяется положение каждого узла новой сетки относительно элементов проецируемой сетки. При нахождении узла внутри элемента значение в нем вычисляется с помощью функций формы соответствующего элемента. Если узел не лежит в элементе, то значение в нем определяется как значение в точке, лежащей в ближайшем элементе. Задание максимального расстояния между узлом и этой точкой позволяет интерполировать значения с заданной точностью. Ниже показана часть командного файла Code_Aster, в которой содержатся операторы, которые считывают сетку в формате *.med и проецируют значения на заданную группу КЭ сетки, содержащую элементы, лежащие на поверхности тела. Полученное поле давлений переводится в тип данных, необходимый для приложения его в качестве нагрузки.
field = LIRE_CHAMP (
TYPE_CHAM = 'NOEU_PRES_R', MAILLAGE = fluid, NOM_MED = 'p', NOM_CMP_MED = '', NOM_CMP='PRES', UNITE = 21 )
proj
PROJ CHAMP
pres
CREA RESU
(
METHODE = 'AUTO', MAILLAGE_1 = fluid, MAILLAGE_2 = mesh,
CHAM_GD = VIS A VIS
field, = F (
TOUT 1
'OUI'
GROUP_MA_2 )
shell'
)
(
OPERATION TYPE_RESU NOM CHAM =
'AFFE', 'EVOL_CHAR' 'PRES',
0)
AFFE= _F (CHAM_GD = proj, INST )
load = AFFE_CHAR_MECA (MODELE = model, EVOL_CHAR = pres)
Недостатком описанной методики является исключение из рассмотрения влияние сил инерции, что не позволяет учитывать существенные динамические перемещения, деформации и напряжения, возникающие при переходных режимах движения и колебаниях конструкции.
2.2. Динамическая постановка
В случае, когда необходимо учитывать динамические нагрузки на лопасть, вызванные силами инерции и пульсациями аэродинамических нагрузок, однако влиянием изменения формы лопасти на характер обтекания можно пренебречь, допустимо использовать методику расчета вынужденных колебаний конструкции, блок-схема которой показана на рис. 2. Для использования методики требуется, чтобы частоты изменения аэродинамических нагрузок и собственные частоты колебаний лопасти были сильно отстроены друг от друга. Методика состоит из следующих этапов:
• проведение аэродинамического нестационарного расчета в пакете ОрепРОЛМ с использованием решателя ртркРОЛМ и определение полей давлений и скоростей как функций от времени;
• перенос полей давлений для каждого момента времени с КО на КЭ сетку и определение нестационарных внешних нагрузок в Code_Aster;
• решение в Code_Aster задачи динамики конструкции под действием полученной нагрузки.
Рис. 2. Алгоритм расчета при динамической постановке задачи Fig. 2. Simulation algorithm at dynamic problem definition
Ниже показан пример скрипта, исполняемого в препостроцессоре Salome-Meca, автоматически преобразующего результаты расчета поля давлений из формата OpenFOAM в формат *.med для некоторого шага (n=5) по времени.
n=5
from paraview import simple as sp
reader = sp.OpenDataFile('/my_directory/field.foam')
sp.SetActiveSource(reader)
sp.Show(reader)
view = sp.GetActiveView()
tsteps = reader.TimestepValues
view.ViewTime = tsteps[n]
sp.Render()
writer = sp.CreateWriter('/my_directory/field.med', reader) writer.UpdatePipeline()
Данная методика не позволяет моделировать режимы автоколебаний.
2.3. Квазистатическая постановка с учётом влияния деформации тела на поток
В случае, когда можно пренебречь динамикой движения лопасти, считать режим обтекания лопасти безотрывным и стационарным, однако требуется учесть влияние упругих перемещений лопасти на условия ее обтекания необходимо определить новое равновесное положение обтекаемой деформированной конструкции в потоке. Для этого применяется более сложная методика расчета, в которой используется итерационный цикл, показанный на рис. 3. Одна итерация поиска нового положения равновесия состоит из следующих этапов:
• проведение стационарного аэродинамического расчета методом контрольного объема в пакете OpenFOAM с использованием решателя simpleFOAM. В результате определяются поля скоростей и давлений;
• перенос результатов аэродинамического расчета из пакета OpenFOAM в Code_Aster;
• определение НДС конструкции методом конечных элементов в Code_Aster;
• передача поля перемещений поверхности конструкции из Code_Aster в OpenFOAM для перестроения или деформации сетки.
Рис.3. Алгоритм расчета при квазистатической постановке задачи с учетом влиянии
деформации тела на поток Fig. 3. Simulation algorithm at quasistatic problem definition taking into account body deformation influence on a stream
Для корректности расчета необходимо, чтобы на каждой последующей итерации поле давлений проецировалось на узлы деформированной от предыдущей нагрузки КЭ сетки. С этой целью на каждой итерации перед проведением процедуры проецирования поля давлений КЭ сетка сдвигалась согласно, рассчитанным перемещениям, после проецирования сетка возвращалась в исходное положение и проводился статический расчет. Для смещения сетки согласно рассчитанным перемещениям могут быть использованы операторы Code_Aster, приведенные ниже.
depl = CREA_CHAMP (
TYPE_CHAM = 'NOEU_DEPL_R', OPERATION = 'EXTR', RESULTAT = stat, NOM_CHAM = 'DEPL', NUME_ORDRE=1 )
mesh = MODI_MAILLAGE (
reuse = mesh, MAILLAGE = mesh,
DEFORME = _F (OPTION = 'TRAN', DEPL = depl) )
3. Применение методик для расчета лопасти ВЭУ 3.1. Исследуемая модель
В качестве исследуемой модели (рис. 4) выбрана лопасть для ветроустановки мощностью 5 МВт [7, 8]. Длина лопасти равна 61,5 метра. Лопасть имеет переменное сечение, состоящее из 6 различных профилей. Площадь поперечного сечения уменьшается от корневого сечения к вершине с максимальным соотношением равным 74,2. Так же лопасть имеет начальную закрутку с плавным увеличением углов установки сечений от вершины к
корню. Максимальный угол закрутки составляет ~ 13,3 градусов. Конструкция лопасти представляет собой оболочку, подкрепленную по всей длине ребрами жесткости.
На практике для изготовления лопастей ВЭУ используются композиционные материалы, к ним относятся стеклоуглеткани, комбинированные армирующие материалы, сотовые плиты и т.п., расчет прочности конструкций из таких материалов является достаточно сложной задачей и является одним из направлений дальнейшего исследования.
Поскольку в данной работе основное внимание уделялось методике объединения пакетов для решения задачи аэроупругости, в расчетах рассматривалась упрощенная модель конструкции, в которой был использован изотропный материал со следующими механическими характеристиками: модуль упругости £=0,68*10п Па, плотность p=2700 кг/м3, коэффициент Пуассона п=0,3. Толщина стенок составила d=1,5 мм.
DU99-W-405
Рис. 4. Исследуемая лопасть ВЭУ Fig. 4. Research blade of the wind turbine
3.2. Расчёт лопасти при квазистатической постановке
Расчет обтекания лопасти в пакете OpenFOAM проводился на неструктурированной гексаэдральной сетке из 400 тысяч элементов с 4 уровнями сгущения сетки.
В качестве граничных условий на лопасти заданы условия жесткой
непроницаемой стенки (для давления - zeroGradient, для скорости - noSlip), а
для боковых граней расчетной области заданы условия свободного потока
(для давления - freestreamPressure, для скорости - freestream со значением 12
м/c). Число Рейнольдса, рассчитанное по характерному размеру тела D
(наибольшая из хорд профилей) составило:
p*lUl*D IUI*D 12 м/с* 5.5 м
ße = ^-= —— = TTr-r-Lm-^rr = 4.37 * 106
ц v 15.1 * 10-6 м2/с
Вычисления проводятся с помощью стационарного решателя simpleFOAM для
несжимаемой вязкой среды. Параметры среды соответствуют параметрам
воздуха при 20 0С (плотность 1,204 кг/м3 и кинематическая вязкость 1,51 *10-5 м2/с). Для моделирования турбулентного течения использовалась модель турбулентности Spalart-AUmaras. Турбулентная вязкость в свободном потоке равняется 0,14 м2/с, на поверхности лопасти 0. Все физические величины в расчетной области определялись в центре расчетной ячейки. Аппроксимация слагаемых в исходных уравнениях была выполнена со вторым порядком точности по времени и пространству. Уравнения для связи скорости и давления решались с помощью итерационного алгоритма SIMPLE [9]. Полученные поля давлений и скоростей среды представлены на рис. 5 и 6. Из рисунков видно, что, когда поток сталкивается с лопастью, он тормозится и изменяет направление движения, обтекая её. При этом около одной поверхности лопасти возникает область с повышенным давлением воздуха, а около другой поверхности - с пониженным. Величина разницы давлений составляет ~ 190 Па. Из-за разности давлений на лопасть начинает действовать аэродинамическая сила. График итерационной сходимости расчета представлен на рис.7, график сеточной сходимости проекций аэродинамической силы, действующей на лопасть представлен на рис. 8.
—
pressure. Pa
-100.0 -50.0 0.00 50.0 100.0
^^^^^^ 1111111111111111111111
Рис. 5. Полученное поле давлений при стационарном расчете Сечения (сверху вниз) x = 12 м, х = 24 м, х = 36 м, х = 48 м Fig. 5. Pressure field at steady-state solution Sections (from top to down) x = 12 m, x = 24 m, x = 36 m, x = 48 m
velocity, m/s
Рис. 6. Полученное поле скоростей при стационарном расчете Fig. 6. Velocity field at steady state-solution
Распределение значения давления на гранях контрольно-объемной сетки показаны на рис. 9.
Для определения НДС лопасти использовалась конечно-элементная модель, состоящая из треугольных оболочечных элементов первого порядка. Размер сетки составил 7714 элементов (из них 1774 приходятся на ребра жесткости). Внешний вид данной сетки показан на рис. 10.
200 400 600
1000 1200 1400 1600 -4662 -4507
Итерация / Iteration
Рис. 7. Сходимости невязки давления, проекций скорости и турбулентной
вязкости Fig. 7. Pressure, speed and turbulent viscosity residuals
Колнчесво ячеек, тыс / Number os cells, ths
—Сила по X —'—Сила no У —•— Сила no Z
Рис. 8. Графики сходимости по сетке значений проекций аэродинамической силы, действующей на лопасть Fig. 8. Aerodynamic force mesh convergence
Рис. 10. Конечно-элементная сетка лопасти Fig. 10. Finite element mesh of the blade Результат применения операций проецирования давления показан на рис.11. Следует отметить, что в данной модельной задаче рассматривалась только аэродинамическая нагрузка, а силы другой природы (центробежные, собственный вес лопасти) не учитывались.
Рис. 11. Поле давлений, интерполированное на КЭ сетку Fig. 11. Pressure field interpolated to the FE mesh
Полученное поле давлений было использовано в качестве внешней нагрузки при определении НДС лопасти, закрепленной по левому торцу. Результаты расчета значений перемещений узлов и напряжений показаны на рис. 12, 13.
Рис. 12. Эпюра перемещений Fig. 12. Plot of the displacement
По результатам расчета можно сделать вывод, что при заданной скорости потока 12 м/с максимальные перемещения на конце лопасти составляют около 0.49 метра, а в конструкции лопасти из изотропного материала не возникает напряжений, приводящих к разрушению.
3.3. Расчёт лопасти в динамической постановке
При решении нестационарной задачи обтекания в OpenFOAM вычисления проводились с помощью нестационарного решателя pimpleFOAM для несжимаемой вязкой среды. В качестве начальных полей давления, скорости и турбулентной вязкости использовались результаты стационарного расчета. Шаг по времени выбирался автоматически из условия для значения числа
Куранта Comax < 0.9. Характерная величина шага по времени - 0.005 секунды. Время расчета переходного режима длительностью 10 секунд составило 2 часа на 1 ядре (Intel(R) Xeon(R) CPU X5670, 2.93GHz). Для проведения расчетов был использован открытый облачный сервис UNIHUB от ИСП РАН [10, 11]. В качестве примера для двух точек, лежащих вблизи поверхности лопасти, на рис. 14 приведены графики зависимости давлений от времени.
Рис. 14. Графики давлений для двух точек (А и В) вблизи поверхности лопасти Fig. 14. Pressure plots for two points (A and B) near the surface of the blade
Полученные в результате аэродинамического расчета поля давлений для каждого момента времени были спроецированы на КЭ сетку, и для каждого момента времени был составлен вектор внешней нагрузки. Динамический расчет переходного режима колебаний лопасти был проведен для той же модели, что и статический. При этом было введено демпфирование по Рэлею с коэффициентами а=Р=0.05. В качестве контролируемого параметра была выбрана У-компонента перемещений точки, лежащей на свободном конце лопасти. Для этой точки был построен график перемещений на интервале времени 0...6 секунды. Шаг интегрирования был выбран 0.1 секунды. На рис. 15 показано сравнение переходных режимов для двух случаев. В первом случае (синий) внешняя нагрузка была нестационарной (менялась на каждом шаге интегрирования), а во втором случае (красный) нагрузка была постоянной, внезапно приложенной в начальный момент времени.
f
1
/ /-Ч У
/
/
У
0 1 2 3 4 5 6
Time, s
Рис. 15. Перемещения точки конца лопасти Fig. 15. Displacement of the tip of the blade
Как видно из графика на рис. 15, в рассмотренном случае учет нестационарности аэродинамической нагрузки поля давлений практически не оказывает значительного влияния на параметры переходного процесса. Объяснить это можно тем, что собственные частоты лопасти (низшая собственная частота 0.56 Гц) значительно ниже частот изменения внешней нагрузки (3.0...4.0 Гц). Максимальное перемещение конца лопасти в переходном режиме составило около 0.58 метра (что на 18% больше статического перемещения в предыдущем расчете), в то время как перемещение в новом статическом положении равновесия, к которому стремиться переходный режим составляет около 0.32 метра, что на 30% меньше чем перемещение в расчете по первой методике.
В дальнейшем расчет может быть проведен с помощью вихреразрешающего моделирования, с использованием метода крупных вихрей, для получения мгновенных составляющих значений скорости и давления [12].
3.4. Расчёт лопасти при квазистатической постановке с учётом влияния деформации конструкции на поток
Для той же самой модели согласно блок-схеме алгоритма, показанной на рис. 3, было проведено 5 итераций (внутренних циклов). Эпюры перемещений лопасти для некоторых итераций даны на рис. 16. График перемещений свободного конца лопасти показан на рис. 17. Видно, что учет влияния деформаций лопасти существенно влияет на параметры ее обтекания, поскольку в данном случае максимальное перемещение составило около 0.53 метра, что на 10% больше, чем в квазистатическом расчете по первой методике.
Рис. 16. Поле перемещения лопасти на каждой итерации Fig. 16. Blade displacement on each iteration
0 1 2 3 4 5
Number of Iteration
Рис. 17. Перемещения точки конца лопасти Fig. 17. Displacement of the tip of the blade
4. Выводы
Работа позволяет сделать вывод о том, что пакеты свободного программного обеспечения OpenFOAM и Code_Aster в сочетании с препостпроцессором Salome-Meca имеют все необходимые средства для создания на их базе единого программного комплекса по решению задач статической аэроупругости и исследованию малых вынужденных колебаний лопасти ВЭУ в потоке вязкой среды.
Предложенные и протестированные на модельной задаче три простых методики связи двух пакетов, основанные на переносе результатов между расчетными сетками, показывают устойчивость счета и адекватность получаемых результатов. Однако величины прогибов, полученные по
различным методикам, существенно различаются, что требует проведения дальнейшей верификации.
В дальнейшем представленные методики могут быть использованы для решения более сложных задач верификации методик по известным экспериментальным данным. Так же полученные результаты являются базой для реализации метода решения полностью связанных задач аэроупругости на базе свободного программного обеспечения.
Благодарности и ссылки на Гранты
Работа выполнена при финансовой поддержке РФФИ (грант № 17-07-01391).
Список литературы
[1]. Бисплингхофф Р.Л., Эшли Х., Халфмэн Р.Л. Аэроупругость. М. Изд-во Иностр. лит. 1958, 800 с.
[2]. Горшков А.Г., Морозов В.И., Пономарев А.Т., Шклярчук Ф.Н. Аэрогидроупругость конструкций. М.: Физматлит, 2000, 592 стр.
[3]. Tukovic Z., Jasak H., Updated Lagrangian finite volume solver for large deformation dynamic response of elastic body. Transaction of FAMENA XXX (1), 2007, pp. 599608.
[4]. Kotsur O., Scheglov G., Leyland P. Verification of modelling of fluid structure interaction fFSIj problems based on experimental research of bluff body oscillations in fluids. In Proceedings: 29th Congress of the International Council of the Aeronautical Sciences, 2014.
[5]. Sekutkovski B., Kostic I., Simonovic A., Cardiff P., Jazarevic V. Three-dimensional fluid-structure interaction simulation with a hybrid RANS-LES turbulence model for applications in transonic flow domain. Aerospace Science and Technology, vol. 49, 2016, pp. 1-16.
[6]. Benra F.-K., Dohmen H. J., Pei J., Schuster S., Wan B. A Comparison of One-Way and Two-Way Coupling Methods for Numerical Analysis of Fluid-Structure Interactions. Journal of Applied Mathematics. Article ID 853560, 2011, p. 16, doi:10.1155/2011/853560
[7]. Resor B. R. Definition of a 5MW/61.5m Wind Turbine Blade Reference Model. SANDIA REPORT SAND2013-2569 Unlimited Release, Printed April 2013, pp.1-53.
[8]. Jonkman J., Butterfield S., Musial W., Scott G. Definition of a 5-MW Reference Wind Turbine for Offshore System Development, National Renewable Energy Laboratory. US, Colorado, Technical Report NREL/TP-500-38060, February 2009, p. 75.
[9]. Weller H.G., Tabor G., Jasak H., Fureby C. A tensorial approach to computational continuum mechanics using object oriented techniques, Computers in Physics, vol.12, № 6, 1998, pp. 620-631.
[10]. Крапошин М.В., Самоваров О.И., Стрижак С.В. Особенности реализации Web-лаборатории механики сплошной среды на базе технологической платформы программы "Университетский кластер". Труды международной суперкомпьютерной конференции, 2011.
[11]. UniHUB: сайт. Режим доступа: http://www.unicluster.ru/unihub.html (дата обращения: 13.11.2017).
[12]. Strijhak S., Redondo J.M., Tellez J. Multifractal analysis of a wake for a single wind turbine. Topical Problems of Fluid Mechanics 2017: Proceedings, 2017. pp. 275-284.
The method of solving aeroelasticity problems for wind blade using open source software
1 2 P.S. Lukashin <ski1l@mail.ru> 1 2 V.G. Melnikova <vg-melnikova@yandex.ru> 2 S.V. Strijhak <strijhak@yandex.ru> 1 G.A. Shcheglov <shcheglov_ga@bmstu.ru> 1 Bauman Moscow State Technical University, 5/1, 2-nd Baumanskaya st., Moscow. 105005, Russia 2 Ivannikov Institute for System Programming of the RAS, 25, Alexander Solzhenitsyn st., Moscow, 109004, Russia
Abstract. Due to the development of Wind Energy and construction of new wind farms in Russian Federation there is a need for the solution of application-oriented problems and development of effective methods for calculation of wind turbine's elements. One of the directions for computational continuous mechanics is connected with problems in aeroelasticity (fluid-structure interaction). The possibility of solving one of the problem in aeroelasticity using a complex program approach on the basis of open source software OpenFOAM and Code_Aster is shown in this article. On the example of the blade for wind turbine, 61.5 meters long, the techniques of solving problem for a static and dynamic aeroelasticity in which calculation of flow of the blade with a subsonic air flow is done in OpenFOAM library (solvers simpleFOAM and pimpleFOAM) are considered. The calculation of the intense deformed status of the blade is done in Code_Aster code. The flowcharts for three different approaches for solving problems of aeroelasticity, examples of scripts and command files for data transfer between two codes in the course of calculation are provided in article. The control-volume mesh consisting their hexahedral elements, the total number is about 400000 elements, for calculation of flow around the blade is constructed in OpenFOAM library, the finite-element mesh consisting of triangular shell elements of first order, the total number is 7714, for calculation of the intense deformed status is constructed in Salome-Meca code. The results of calculation are provided in the form of fields for pressure and velocities; graphics for residuals of pressure, velocity, turbulent viscosity; projections of aerodynamic force from time; diagrams of displacement and stress; the values of pressure for two points for the surfaces and displacement of the tip of the blade from time. The calculations are run using resources of UniHUB web-laboratory ISPRAS.
Keywords: aeroelasticity; wind blade; solver, Code_Aster; OpenFOAM; finite-volume method, finite-element method, mesh, turbulence model, static, dynamic, frequency, displacement, stress, calculations, web-laboratory.
DOI: 10.15514/ISPRAS-2017-29(6)-16
For citation: Lukashin P.S., Melnikova V.G., Strijhak S.V., Shcheglov G.A. The method of solving aeroelasticity problems for wind blade using open source software. Trudy ISP RAN/Proc. ISP RAS, vol. 29, issue 6, 2017, pp. 253-270 (in Russian). DOI: 10.15514/ISPRAS-2017-29(6)-16
References
[1]. Bisplinghoff R.L., Ashley H., Halfman R.L. Aeroelasticity. Dover Publications, Inc., Mineola, N.Y., 1996, 800 p.
[2]. Gorshkov A.G., Morosov V.I., Ponomarev A.T., Schklyaruk F.N. Aeroelasticity of structures. M.: Fizmatlit, 2000, 592 p. (in Russian).
[3]. Tukovic Z., Jasak H., Updated Lagrangian finite volume solver for large deformation dynamic response of elastic body. Transaction of FAMENA XXX (1), 2007, pp. 599608.
[4]. Kotsur O., Scheglov G., Leyland P. Verification of modelling of fluid structure interaction (FSI) problems based on experimental research of bluff body oscillations in fluids. In Proceedings: 29th Congress of the International Council of the Aeronautical Sciences, 2014.
[5]. Sekutkovski B., Kostic I., Simonovic A., Cardiff P., Jazarevic V. Three-dimensional fluid-structure interaction simulation with a hybrid RANS-LES turbulence model for applications in transonic flow domain. Aerospace Science and Technology, vol. 49, 2016, pp. 1-16.
[6]. Benra F.-K., Dohmen H. J., Pei J., Schuster S., Wan B. A Comparison of One-Way and Two-Way Coupling Methods for Numerical Analysis of Fluid-Structure Interactions. Journal of Applied Mathematics. Article ID 853560, 2011, p. 16, doi:10.1155/2011/853560
[7]. Resor B. R. Definition of a 5MW/61.5m Wind Turbine Blade Reference Model. SANDIA REPORT SAND2013-2569 Unlimited Release, Printed April 2013, pp.1-53.
[8]. Jonkman J., Butterfield S., Musial W., Scott G. Definition of a 5-MW Reference Wind Turbine for Offshore System Development, National Renewable Energy Laboratory. US, Colorado, Technical Report NREL/TP-500-38060, February 2009, p. 75.
[9]. Weller H.G., Tabor G., Jasak H., Fureby C. A tensorial approach to computational continuum mechanics using object oriented techniques, Computers in Physics, vol.12, № 6, 1998, pp. 620-631.
[10]. Kraposhin M.V., Samovarov O.I., Strijhak S.V. The features of design Web-laboratory of computational continuum mechanics on the basis of the technological platform in scope of the program "University Cluster". The Proceedings of the International Supercomputer Conference, 2011 (in Russian).
[11]. UniHUB. Available at: http://www.unicluster.ru/unihub.html (Accessed 13 November 2017).
[12]. Strijhak S., Redondo J.M., Tellez J. Multifractal analysis of a wake for a single wind turbine. Topical Problems of Fluid Mechanics 2017: Proceedings, 2017. pp. 275-284.
