Как декантировать вино. ANSYS CFD на службе у сомелье.

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

Для тех, кто ведет трезвый образ жизни, напоминаю, что декантирование (декантация) — это процесс отделения вина от осадка и насыщение его кислородом. При этом сомелье совершает руками серию магических движений. А сама процедура может занимать 10-15 минут и более. Чаще всего декантируют красные вина и реже белые. Отметим, что споры о полезности и объективности результатов этого процесса не утихают до сих пор. Поэтому работа японских инженеров заслуживает внимания. Они провели исследование бокалов различной формы на предмет интенсификации (или сокращения) процесса аэрации вина при контакте с кислородом.

2016-11-02_10-13-15

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

bokaly

В 17-м релизе ANSYS существенно упростилась процедура постановки 6DOF-задачи. Раньше для задания движения тела с 6-ю степенями свободы  необходимо было создавать и компилировать специальные UDF-макросы.  Но в Fluent 17.0 появилась специальная панель, которая сделала работу инженера-исследователя намного легче. При этом сохранилась возможность определять закон движения тела как на основе собственных пользовательских функций (UDF), так и с использованием возможностей System Coupling.

6dof

По результатам моделирования в номинации «Лучший аэрационный бокал-2016» победил бокал на длинной утонченной ножке.

rezultatyЕсли вы в совершенстве владеете японским языком, то более подробно про это исследование вы можете почитать по этой ссылке http://news.mynavi.jp/kikaku/2015/12/16/004/. А к капиллярам мы вернемся уже в следующем году.

С уважением, Денис П. Хитрых,
Директор АО «СимуЛабс».
2016-11-21_16-58-13

Мастер-класс: CFD расчет химического реактора #4.

reactorСегодня я опубликую последнюю (крайнюю) часть моего мастер-класса по расчету химического реактора в ANSYS CFX. Ранее мы с вами научились определять многоступенчатые реакции в CFX на основе предопределенных шаблонов; построили геометрическую модель реактора; сгенерировали расчетную сетку и пр. И нам осталось только определить граничные условия и настройки решателя. Чем мы сейчас и займемся.

На входе мы задаем массовый расход, параметры турбулентности и температуру реагентов.

2016-12-26_17-04-54

2016-12-26_17-05-42

Турбулентность мы определим через относительную интенсивность, равную 0.05, что эквивалентно Tu = 5%. Массовый расход = 0.6 кг/с и статическая температура = 333 К.

Кроме того, необходимо задать компонентный состав на входе в соответствии с таблицей, показанной ниже.

2016-12-26_17-26-55

Если просуммировать массовые доли всех компонентов на входе, то мы получим величину, равную 0.9868 (вместо 1). Пугаться не стоит: «недостающиеся» доли  (1-0.9868 = 0.0132) приходятся на «замороженный» (constraint) азот, для которого не решается уравнение переноса (см. рис. ниже).

Для определения «замороженного» компонента вам необходимо отредактировать свойства расчетного домена. Откройте закладку Fluid Models и укажите для азота N2 опцию Constraint. Для всех других компонент по-умолчанию будет установлена опция Transport Equation.

transport

Не выходя из режима редактирования свойств расчетного домена, задайте величину опорного давления (Reference Pressure) равной 0.35 МПа. В закладке Fluid Models для Heat Transfer укажите Thermal Energy; Turbulence = SST и Combustion = FRC (модель, описывающая скорости химических реакций — «антагонист» модели «быстрой» химии). В Material выбираем ранее созданный материал Chloroform.

На выходе определяем статическое давление, равное 0.0 Па. И, наконец, на стенках определяем граничное условие первого рода T = const = 800 К.

Для начальной инициализации оставляем все настройки по-умолчанию. Задаем максимальное количество итераций = 400. Шаг по времени = 0.1 сек. И критерий сходимости RMS = 1e-04.

На этом все. Запускаем задачу на решение.

С уважением, Денис П. Хитрых,
Директор АО «СимуЛабс».
2016-11-21_16-58-13

Мастер-класс: CFD расчет химического реактора #3.

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

Геометрия расчетной области в виде цилиндра показана на рисунке ниже. Все размеры даны в относительных величинах, привязанных к характерному размеру d = 10 мм. В статье, которую мы взяли в качестве основы для решения нашей задачи, выполнено моделирование для реактора с предсмешением, в котором два потока с разным химическим составом подводятся через кольцевой канал и круглую трубу диаметром 12 мм, соответственно. А вся геометрическая модель образмеривается относительно диаметра этой внутренней трубки.

В нашем случае никакого предсмешения нет, и поток поступает в зону реакции через обычный круглый канал диаметром 10 мм. Еще одно отличие нашей геометрии от референсной заключается в длине реактора. Мы искусственно увеличили ее почти в 2 раза. Все другие размеры практически совпадают.

geometry

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

Какие сложности могут у вас возникнуть при генерации расчетной сетки?  Рассмотрим вариант задачи с циклосимметрией. На рисунке ниже я показал основные проблемы, которые могут возникнуть при разбиении сетки в режиме «по-умолчанию».

mesh

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

Вторая проблема — это наличие призматических элементов на оси реактора по всей его длине. Избавится от  призм можно двумя способами. В первом случае необходимо с помощью виртуальных узлов, ячеек, ребер и пр. (Virtual Topology) построить L-топологию на входном и выходном каналах, соответственно . Другой вариант — вручную декомпозировать расчетную область на sweepable-объемы, как показано на рисунке ниже.

sweepable

Когда все проблемы решены, можно смело нажимать кнопку [Generate].

С уважением, Денис Хитрых,
АО «СимуЛабс».
2016-11-21_16-58-13

Мастер-класс: CFD расчет химического реактора #2.

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

Продолжаем заниматься расчетами химических реакторов в ANSYS CFX. Прошлый наш разговор мы закончили описанием инструментов  ANSYS CFX для задания элементарных реакций. Напомню, что для вызова панели [Reaction] используется иконка reaction1, которая находится в верхнем контекстном меню.  Для определения реакции (Insert Reaction и далее вбиваете имя реакции) вы последовательно проходите все закладки этой панели (слева-направо): Basic Settings→ Reactants→ Products→ Reaction Rates.

Чтобы упростить вам жизнь, я заранее подготовил CCL-файл с библиотекой всех реакций. Для его активации вы должны выполнить команду File→ Import→ CCL… Файл с реакциями вы можете скачать по этой ссылке

.

В закладке Reactants мы указываем исходные реагенты, в закладке Products — конечные реагенты (продукты реакции). А в закладке Reaction Rates мы определяем кинетику соответствующей реакции. Все реагенты, которые участвуют у нас в реакциях (а это начальные, промежуточные и конечные продукты), мы должны заранее определить с помощью инструмента materialMaterial . После этого все действующие реагенты нам станут доступны через список Materials List.

С помощью инструмента Material мы определяем теплофизические свойства реагентов (как и в обычной газодинамической задаче). Это достаточно трудоемкий процесс, так как нам необходимо создать 16 новых «материалов»: Cl2, Cl, CH3Cl, CH2Cl, HCl, CH2Cl2, CHCl2, CHCl2CH2Cl, CCl3, CCl3CCl3, CHCl2CHCl2, CHCl3, CCl3CH2Cl, CCl3CHCl2, CCl4, N2). Пример задания свойств хлора показан на рисунке ниже.

clГотовую библиотеку материалов для проекта, сохраненную в формате CCL, вы можете скачать здесь

.

«Подключаете» к проекту вы ее также, как и библиотеку реакций.

На сегодня это все. Я думаю, что в следующий раз мы займемся построением расчетной CAD-модели и и поговорим о расчетной сетке для данной задачи.

С уважением, Денис Хитрых,
АО «СимуЛабс».
2016-11-21_16-58-13

Мастер-класс: CFD расчет химического реактора

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

Упрощенный (редуцированный) химический механизм превращения хлорметана в хлороформ включает 16 реагентов (исходных и конечных) и 21 реакцию, которые показаны в таблице ниже. В этой же таблице в двух крайних правых столбцах представлены значения для энергии активации E и предэкспоненциального фактора A, необходимые для задания температурной зависимости констант скоростей элементарных реакций в препроцессоре ANSYS CFX (или в другом газодинамическом пакете) при описании цепного механизма реакций.

reaction-mechanism_21_reaction

Более детально химический процесс образования хлороформа описывается 152 реакциями, в которых участвуют 38 реагентов (химических соединений). Описание этого механизма вы можете найти в работе Jimmy J. Shah  и Rodney O. Fox из Iowa State University: «Computational Fluid Dynamics Simulation of Chemical Reactors…» (1999).

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

  1. Что изучает химическая кинетика? Химическая кинетика изучает закономерности развития химических реакций.
  2. Почему описание даже относительно простого  химического процесса, например, горения водорода включает почти 40 элементарных реакций? Формально горение водорода выражается одной глобальной реакцией: 2H2 + O2 = 2H2O. Т. е. можно предположить, что в результате одного столкновения этих трех молекул (двух молекул водорода и одной молекулы кислорода) образуется две молекулы воды. С позиции теории вероятности такое событие является маловероятным (точнее невозможным). Поэтому реакция горения водорода протекает поэтапно (в несколько стадий) через промежуточное образование атомов водорода, кислорода и радикалов.  Реакции такого типа называются цепными или суммарными. Цепные реакции являются последовательностью большого числа элементарных реакций, выявлением и изучением которых занимаются химики и физики.
  3. Стехиометрическое уравнение (или уравнение реакции). Стехиометрическое уравнение представляет собой краткое выражение материального баланса реакции. Например, уравнение 2H2 + 1O2 = 2H2O означает, что всякий раз, как в процессе реакции затрачиваются две молекулы водорода, одновременно расходуется ровно одна молекула кислорода и образуются две молекулы воды. Коэффициенты перед реагентами называются стехиометрическими коэффициентами. Если реакция состоит из ряда стадий, то получается система из n стехиометрических уравнений.
  4. Закон Аррениуса. Для химических реакций характерна сильная нелинейная зависимость констант скорости k от температуры. Эта температурная зависимость описывается достаточно простой формулой — законом Аррениуса: k = Aexp(-Ea/RT).  Здесь Ea — это энергия активации (она имеет размерность [Дж/моль] и A — предэкспоненциальный фактор (или частотный фактор). Размерность A совпадает с размерностью k. Обратите внимание, что в ANSYS CFX  предэкспоненциальный фактор имеет размерность [time^-1 (mol m^-3)^(1-n)], где n — суммарный порядок реакции.
  5. Быстрая или медленная химия. В теории горения есть критерий подобия — число Дамкелера Dm, который определяет отношение скорости течения химической реакции к скорости физического процесса. Если Dm→0, то время протекания химической реакции намного больше характерного времени физического процесса. Т. е. можно не учитывать химические реакции и считать смесь газов химически инертной (как-бы «замороженной»). Другими словами, реакция развивается намного медленнее (среда не успевает измениться) по сравнению с изменениями гидродинамических параметров. Второй предельный случай противоположен первому (Dm→∞): в каждой точке потока очень быстро устанавливаются такие концентрации реагентов, которые соответствуют равновесному составу. Оценив значение числа Дамкелера, можно в первом приближении выбрать соответствующую модель горения: для «быстрой» химии — это модель EDM, для «медленной» химии — это модель FRC.

Читать далее

Место встречи двух океанов. Разрушаем иллюзии.

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

mixing

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

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

При контакте двух неподвижных взаиморастворимых жидкостей появляется достаточно узкая переходная зона (transition zone), которая постепенно увеличивается в размерах за счет процесса диффузии. При существенных различиях в плотности жидкостей может возникнуть конвективное течение со слабовыраженными «периодическими» структурами, которое часто ошибочно интерпретируют как неустойчивость Рэлея-Тейлора.

mixing2В нашем эксперименте две жидкости (пресная и солёная вода) заливались в общую ванну, разделенную подвижной вертикальной перегородкой на две равные по объему части. Затем перегородку медленно выдвигали и наблюдали процесс «смешения» пресной и соленой воды.

При использовании метода VOF (метод объема жидкости) мы сможем смоделировать динамику свободной поверхности —  межфазной границы раздела двух сред.

В методе VOF для определения положения свободной поверхности решается дополнительное уравнение конвективного переноса объемной доли жидкости в ячейке сетки, которая выступает в качестве т. н. функции «маркера» (и принимает значения в интервале от 0 до 1: ячейка заполнена жидкостью или ячейка пуста).

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

Читать далее

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.