Как «смоделировать» ветрогенератор в ANSYS CFD

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

Поэтому на следующей неделе мы рассмотрим с вами наиболее интересные варианты конструкций ветряных турбин. И научимся моделировать в ANSYS CFD аэродинамику классического осевого «ветряка».

Для начала немного истории. На фотографиях снизу изображен один из первых отечественных ветрогенераторов времен СCCP. Он был спроектирован в ЦАГИ и введен в эксплуатацию в 1931 году. Мощность станции — 100 кВт (при скорости ветра 10 м/с); место установки — Балаклава, Крым. В 1943 году ветряная станция была разрушена. Сейчас на высоте Горной можно найти следы фундамента опор этой станции.
turbina

По современным меркам 100 кВт — это немного, но в начале 30-х годов прошлого века эта цифра заслуживала уважения. Например, немцы и датчане в это же время строили ветроэлектростанции меньшей мощности и с меньшим диаметром лопастей. Примечательно, что мачта Балаклавской станции была построена по проекту Владимира Шухова (автора знаменитых строительных гиперболоидов). А сама станция питала энергией Севастопольский трамвайчик.

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

windturbine

А это турбина, которую мы с вами рассчитаем в ANSYS CFX или ANSYS Fluent. Размер расчетной области здесь показан условно.

До скорой встречи.

С уважением, Денис Хитрых (АО «СимуЛабс»).

Инфографика: Сопло Лаваля

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

Немного о сопле Лаваля
Чтобы разогнать газ от полностью заторможенного до сверхзвукового необходим специальный канал, который сначала сужается, а затем расширяется. В самом узком сечении этого канала, число Маха должно равняться единице. Такой канал получил название сопла Лаваля. Особенности режимов течения для сопла Лаваля можно установить с помощью упрощенных диаграмм. Для этого необходимо знать соотношение параметров торможения на входе в сопло и параметров среды, в которую истекает струя на выходе из сопла. Зная параметры торможения на входе и закон профилирования сопла вдоль его оси можно рассчитать газодинамические параметры потока по всей длине сопла.

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

Мы надеемся, что объем информации, заложенный в нашу инфографику, поможет вам сделать первые шаги на пути освоения технологий вычислительной гидродинамики (CFD). И первым вашим смоделированным техническим устройством будет именно сопло Лаваля.

Что такое сопло Лаваля

На сайте размещена черновая версия инфографики. Все замеченные опечатки устранены в печатной версии плаката. Первые 200 копий инфографики будут бесплатно розданы всем участникам пользовательской конференции ACUM-2016 (Москва, 25-27 октября). Вы также можете подъехать за плакатом к нам в офис, который располагается на ул. Суздальская, дом. 46.

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

С уважением, Денис Хитрых (АО «СимуЛабс»).

1 апреля — День ANSYS Fluent на cfd-blog.ru

ВелосипедСегодня мы целый день посвятим расчетному CFD-пакетуANSYS Fluent. И рассмотрим с вами около двух десятков вопросов, которые поступили в адрес службы технической поддержки компаний ANSYS и КАДФЕМ.

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

fluentСамый первый вопрос связан с инсталляцией ANSYS Fluent версии 17.0. На моей рабочей станции установлена операционная система Windows 7 и при инсталляции 17-й версии ANSYS у меня (и, вероятно, у многих других пользователей) возникли проблемы при работе c GUI Fluent. Решение данной проблемы заключалось в установке пакета обновлений для Window 7 — KB2533623, который предназначен для повышения безопасности работы со внешними библиотеками. Скачать этот пакет можно с официального сайта Microsoft по следующей ссылке https://support.microsoft.com/ru-ru/kb/2533623.

fluentСледующий вопрос связан с оценкой сходимости нестационарной задачи в ANSYS Fluent на каждом временном шаге по заданному критерию сходимости (например, по изменению величины какого-либо интегрального параметра). Для этого можно использовать команду stptmstp (=SToP TiMe STeP) языка Scheme. Более детальное описание этой команды можно получить набрав в консоли (Console) следующую последовательность символов: (stptmstp-docu-alt).

Все значения переменных, которые используются в команде stptmstp по умолчанию, можно отредактировать. Например, набрав команду вида (set! stptmstp-n 7), вы измените количество итераций, за которые собирается статистика по заданной переменной, с 3-х (это значение по умолчанию) до 7-и.

Или, например, с помощью команды set! stptmstp-maxrelchng вы можете изменить значение  относительной разницы между собственными значениями заданной переменной, вычисленными в последовательных итерациях. По умолчанию, ее величина соответствует 1.e-03. Таким образом, для того, чтобы изменить значение этой переменной с 1.e-03 на 1.e-5, вам следует набрать команду set! stptmstp-maxrelchng 1.e-05.

Далее мы приведем небольшой пример использования команды stptmstp в интерфейсе Fluent. Откройте последовательно закладки Solution → Calculation Activitie s→ Execute commands. Сначала введите первую команду command-1: (stptmstp-resetvalues). И выберите опцию Time Step. Эта команда «сбросит» анализ предыдущего шага по времени. После этого введите команду  (stptmstp-chckcnvrg ‘ ( «/report/surface-integrals/integral wall () heat-flux n»)) и выберите опцию Iteration. С помощью этой команды будет выполнен анализ последних трех итераций (3 итерации — это значение по умолчанию в ANSYS Fluent, как мы уже говорили выше), и будет выполнена проверка насколько изменился тепловой поток на заданной поверхности за эти три последних итерации. При этом относительное изменение ограничено значением 1e-03 (по умолчанию).

execute

fluentТретий вопрос регулярно обсуждается на различных CFD-форумах, его также часто адресуют в нашу службу технической поддержки. Но сегодня, чтобы дать на него полноценный ответ, мы воспользуемся популярным ресурсом Eureka. Задача формулируется следующим образом: как с помощью UDF определить табулированные свойства пользовательского материала. Соответствующий шаблон UDF приведен ниже (в раскрывающемся списке). Вы можете легко доработать его для решения вашей конкретной задачи.

Код UDF

#include "udf.h"

int NPts_mu = 10; /* number of entries in viscosity table */

/* Locate the place in the interpolation vector using bisection algorithm*/
int locate(float xx[], int n, float x)
{
int j = 0;
int jm = 0;
int jl = 0;
int ju = n+1;
int ascnd = 0;

ascnd = (xx[n] >= xx[1]);
while (ju-jl > 1)
{
jm = (ju+jl) >> 1;
if (x >= xx[jm] == ascnd)
jl = jm;
else
ju = jm;
}

if (x == xx[1])
j = 1;
else if (x == xx[n])
j = n-1;
else
j = jl;

return j;
}

float *temp_vec_mu,*visc_vec;
#define FAST_LOOKUP TRUE /* use bisection algorithm for interpolation */
#define TABLE_MESSAGES TRUE
#define DISPLAY_TABLES TRUE

/* Obtaine mu given temperature */
float get_mu_from_T(float xdata)
{
int i = 0;
int j = 0;
float xL,xU,ydata;

#if FAST_LOOKUP
j = locate(temp_vec_mu,NPts_mu,xdata);
xL = temp_vec_mu[j];
xU = temp_vec_mu[j+1];
ydata = visc_vec[j] + (xdata-xL)/(xU-xL)*( visc_vec[j+1] - visc_vec[j] );
#else
for ( i=1; i=xL)&&(xdata<=xU) ) { ydata = visc_vec[i] + (xdata-xL)/(xU-xL)*( visc_vec[i+1] - visc_vec[i] ); break; } } #endif if ( xdata>temp_vec_mu[NPts_mu] )
{
#if TABLE_MESSAGES
Message("n temperature is above the bound of visc-array n");
#endif
ydata = visc_vec[NPts_mu];
}
if ( xdata<temp_vec_mu[1] )
{
#if TABLE_MESSAGES
Message("n temperature is below the bound of visc-array n");
#endif
ydata = visc_vec[1];
}

return ydata;
}

/* Read in the data file containing viscosity as a function of temperature */
DEFINE_ON_DEMAND(read_viscosity)
{
int i = 0;
float xdata,ydata;
FILE* fp;

fp = fopen("viscosity.dat","r");
if ( fp!=NULL )
{
#if !RP_NODE
Message(" n");
Message("Reading file viscosity.dat n");
Message(" n");
#endif
}
else
{
#if !RP_NODE
Message(" n");
Message("Error opening file n");
Message(" n");
#endif
}

temp_vec_mu = (float *) malloc(NPts_mu*sizeof(float));
visc_vec = (float *) malloc(NPts_mu*sizeof(float));

if ( (temp_vec_mu==NULL)||(visc_vec==NULL) )
{
#if !RP_NODE
Message("Memory allocation error n");
#endif
}

for ( i=1;i<=NPts_mu;i++ )
{
fscanf(fp,"%f %e n",&xdata,&ydata);
temp_vec_mu[i] = xdata;
visc_vec[i] = ydata;
#if DISPLAY_TABLES
#if !RP_NODE
Message("%.1f %e n",temp_vec_mu[i],visc_vec[i]);
#endif
}
#else
}
#if !RP_NODE
Message(" n");
#endif
#endif

fclose(fp);
}

/* Interpolate viscosity from look-up table and assign to cell value */
DEFINE_PROPERTY(viscosity,c,t)
{
float mu_lam;
float temp = C_T(c,t);

/* interpolate mu_lam */
mu_lam = get_mu_from_T(temp);

return mu_lam;
}


fluentЧетвертый вопрос будет касаться начальной инициализации задачи со свободной поверхностью в ANSYS Fluent. Самый простой вариант — определить свободную поверхность как плоскую поверхность. Но, что делать, если геометрия свободной поверхности имеет более сложную форму?

Есть несколько вариантов решения данной проблемы. Вы можете использовать специальный «генератор» волн, встроенный в Fluent — Open Channel Wave Model (волны в открытом канале). При инициализации границы свободной поверхности (Free Surface Level) вместо опции Flat укажите опцию Wavy (см. Open channel Init. Method).

wavy

Если вы проводите инициализацию с помощью TUI, то после выполнения команды solve → initialize → open-channel-auto-init, вам необходимо выбрать опцию wavy free surface initialization.

Для инициализации границы свободной поверхности вы также можете использовать Custom Field Functions. Это своеобразный аналог CEL-выражений пакета ANSYS CFX.

VOF

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

Более подробно о двух последних вариантах я расскажу в другой раз.

На сегодня всё. Поздравляю всех с праздником Всемирного Дурака, а в качестве подарка предлагаю вам полюбоваться на обновленный дизайн графического интерфейса ANSYS Fluent 17.0 (см. рис. ниже).

Измененный интерфейс ANSYS Fluent.

Измененный интерфейс ANSYS Fluent.

Традиционный интерфейс ANSYS Fluent.

Традиционный интерфейс ANSYS Fluent.

На следующей неделе мы продолжим отмечать 1 апреля вместе с ANSYS Fluent. Удачного Вам моделирования!

KhitrykhDP

 

С уважением, Денис Хитрых,
Директор SimuLabs4D.

ДР С. П. Королёва — главного ракетостроителя СССР и 20-го века. Гиперзвук, скачок уплотнения и ПС. Моделируем в ANSYS CFX

SPKorolevСегодня мы отмечаем день рождения Сергея Павловича Королёва. Так получилось, что отечественные СМИ в очередной раз проигнорировали это важное событие.

Для тех, кто родился в 90-е и 00-е годы прошлого и текущего столетия, хочу напомнить, что без непосредственного участия этого человека сегодня не работало бы ни спутниковое телевидение, ни системы спутниковой навигации GPS/ГЛОНАСС и многие другие полезные вещи.

Поздравляем коллективы ОАО РКК «Энергия», ФГУП ЦНИИМАШ и других сопричастных организаций с этим замечательным событием. Желаем почаще смотреть на звёздное небо и искать на нём источники своего вдохновения.

Учитывая, что праздник напрямую связан с космосом, задача, которую мы сегодня с  вами рассмотрим, так же будет «космической».

В качестве объекта моделирования мы выбрали модель в виде цилиндра, сопряженного с конусом (см. рис. ниже). Это эталонная модель, которая используется при тестировании моделей турбулентности при различных режимах обтекания исследуемого объекта. В нашей задаче М = 7.05, что соответствует  гиперзвуковому режиму обтекания.

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

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

hyper

Для данной модели есть экспериментальные данные, полученные Kussoy M. I. и Horstman C. C. в 1989 г. (Documentation of Two- and Three- Dimensional Hypersonic Shock Wave / Turbulent Boundary Layer Interaction Flows. NASA Technical Memorandum 101075). Ссылка для скачивания http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19890010729.pdf.

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

Читать далее

Моделирование движения твердого тела в ANSYS CFX c использованием Rigid Body Solver

Сегодня мы кратко рассмотрим задачу моделирования движения твердого тела в ANSYS CFX c использованием Rigid Body Solver. Так как в общем случае твердое тело обладает шестью степенями свободы, то и общая система уравнений движения твердого тела содержит шесть независимых уравнений. Эти уравнения можно представить в виде двух векторных уравнений для скорости изменения импульса и момента импульса тела:

Equation

Здесь P — это полный импульс тела; F — сумма всех сил, действующих на тело со стороны внешних источников; π — момент импульса тела; m — сумма моментов всех сил, действующих на тело.

Mомент импульса π связан с угловой скоpостью посpедством тензоpа инеpции Ι:

moment

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

Уравнения поступательного и вращательного движения твердого тела имеют следующий вид:

equ_1

Здесь FAero — это аэродинамическая сила; FExt — другие внешние силы, действующие на тело; kSpring — линейная жесткость пружины; mAero — аэродинамический момент; kRotate — крутильная жесткость пружины.

На этом мы завершим наш короткий теоретический раздел (более подробно см. документацию ANSYS CFX) и перейдем к практической части.

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

missle

Читать далее

Аэродинамика. Аэроакустика. Часть 3

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

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

Иногда на небе мы можем наблюдать облака странной и оригинальной формы — это т. н. волны Кельвина-Гельмгольца. Это классический пример сдвиговой неустойчивости (неустойчивость тангенциальных разрывов скорости), которая реализуется, когда один слой жидкости «скользит» по другому. Подобная ситуация свойственна многим реальным природным течениям. Неустойчивости в данном случае подвержено положение границы между слоями жидкости, которые движутся с различными скоростями.

kelvin

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

Если говорить о постановка нестационарной задачи, то она имеет много общего с постановкой стационарной задачи. К отличиям можно отнести большой объем расчетных данных, специальные методы для обработки этих данных и специфические настройки CFD-решателя.

Какие вопросы мы должны решить при постановке нестационарной задачи в Fluent?

  1. Как поставить и запустить на расчет нестационарную задачу.
  2. Как выбрать правильный шаг по времени для нестационарной задачи.
  3. Как правильно обработать результаты расчета нестационарной задачи.

tsНачнем с шага по времени. Основное требование — величина шага по времени должна коррелировать с характерным временем моделируемого процесса/явления.

Для стандартной задачи шаг по времени можно оценить из следующего соотношения: delta_t = L/3*V, здесь L — характерная длина, V — характерная скорость.

Для лопаточной машины:  delta_t = N/10*omega, где N — количество лопаток (лопастей), omega — скорость вращения.

Для задачи естественной конвекции: delta_t = L/(g*beta*delta_T*L)^1/2. Здесь g — ускорение свободного падения; beta -температурный коэффициент объемного расширения; delta_T -разница температур (стенка-жидкость).

Наконец, для задачи теплопроводности: delta_t = L^2/(lambda/(rho*Cp)), где lambda и Cp — теплопроводность и теплоемкость среды, соответственно.

Продолжение следует…

Удачного Вам моделирования!

Денис Хитрых, директор R&D Центра SimuLabs4D.

Аэродинамика. Аэроакустика. Часть 2

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

Первое. Что такое средняя длина хорды (СДХ)? В традиционной литературе используется более корректный термин  — «средняя аэродинамическая хорда» (см. рис. ниже). Далее мы будем пользоваться также обозначением Сmac (Mean Aerodynamic Chord ).

хорда

Второе. Я указал требуемое количество элементов в пристеночной области, но ничего не сказал об общей размерности подобных задач, а также геометрии и размерах расчетной области. Если говорить о задаче нестационарной аэродинамики самолета (в конфигурации «фюзеляж+крыло+мотогондола»), то начинать расчет можно с сетки размерностью около 15-20 млн. ячеек. Далее следует увеличивать суммарное количество ячеек с коэффициентом 2,5-3,2. При этом следует учитывать, что с изменением объема сетки, будут изменяться параметры сетки и в пристеночной области. Т. е. с ростом суммарного количества ячеек значение y+ будет уменьшаться. Например, в требованиям к расчетной сетке для 1-го AIAA CFD High Lift Prediction Workshop (2009 г.) указано, что для Re = 4,3e6 и Сref = Cmac =1,00584 м (39,6 дюйма) шаг по высоте первой пристеночной ячейке задается следующей последовательностью: y+ ~ 1,0 -> dy = 0,00020 дюйма;  y+ ~ 2/3 -> dy = 0.00013 дюйма; y+ ~ 4/9 -> dy = 0.00009 дюйма; y+ ~ 8/27 dy = 0.00006 дюйма.

Если говорить о коэффициенте роста элементов, то рекомендуется использовать значения, не превышающие величины 1,25 (для любого типа сетки). На кромки крыла рекомендуется наносить от 6 до 20 ячеек. Шаг сетки по размаху крыла — 0,1% от величины полуразмаха крыла.

Теперь к вопросу о расположении входной границы (например, Farfield в Fluent). В упомянутом выше воркшопе поверхность Farfield располагалась на расстоянии ~100 Cref от планера.

Третье. Возвращаясь к теме выбора значения ШГВ, отмечу, что все приведенные значения относятся к варианту с условно точной сеткой, в которой y+ лежит в диапазоне от 0,4 до 0,5.

Четвертое. Где можно получить дополнительную информацию о  модельмодели турбулентности? Найти описание этой модели можно в стандартной  документации ANSYS или в статье (Transition Modeling for General CFD Applications in Aeronautics. R.B. Langtry and F.R. Menter), которую можно скачать здесь

.

Если говорить о калибровке данной модели, то обычно калибруются две эмпирические поправочные функции:  Fl и Re, соответственно. Первая контролирует длину области перехода. Вторая — это критическое число Рейнольдса, которая определяет начало роста перемежаемости в пограничном слое. Для тех, кто незнаком с этим термином, могу сказать, что под этим термином понимают локальное нарушение однородности турбулентности, когда активные области сосуществуют с условно ламинарными областями.

На этом все. В третьей части нашего мини-семинара мы поговорим о настройках решателей CFX и Fluent, методе LES, подсеточных моделях и пр.

С уважением, Денис Хитрых,

R&D Центр SimuLabs4D.

www.simulabs.ru