Нейросетевое моделирование неизотермических течений неньютоновскихжидкостей Текст научной статьи по специальности «Физика»

УДК 536.24

НЕЙРОСЕТЕВОЕ МОДЕЛИРОВАНИЕ НЕИЗОТЕРМИЧЕСКИХ ТЕЧЕНИЙ НЕНЬЮТОНОВСКИХ ЖИДКОСТЕЙ С.Г. Валюхов, А.В. Кретинин, О.В. Стогней, В.В. Костенко

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

Ключевые слова: гидродинамика, неньютоновские жидкости, нейросетевое моделирование

1. Введение

Явный метод маркеров и ячеек МАС [1] может быть адаптирован для моделирования тепломассопереноса жидкостей, согласующихся со степенным реологическим уравнением [2, 3]. Центральным звеном классического метода маркеров и ячеек является численное решение уравнения Пуассона. При конечно-разностной аппроксимации производных в приграничных расчетных узлах со вторым порядком точности не обеспечивается диагональное преобладание матрицы коэффициентов системы сеточных алгебраических уравнений, что делает в общем случае открытым вопрос применения итерационных процедур ее решения. Использование прямых методов решения систем сеточных линейных уравнений большой размерности ограничено проблемой численной диффузии. Кроме того, псевдонестационарный алгоритм МАС при реализации методом установления итераций по времени плавно перестраивает решение для распределения давления в области течения в зависимости от соответствующего изменения правой части уравнения Пуассона. Эти особенности метода МАС побудили исследовать возможность применения для решения уравнения Пуассона метода взвешенных невязок на базе нейросетевых пробных функций (НМВН) [4-6].

Валюхов Сергей Георгиевич - ВГТУ, д-р техн. наук, профессор, тел. (473) 252-34-52

Кретинин Александр Валентинович - ВГТУ, д-р техн. наук, доцент, e-mail avk-vrn@mail.ru

Стогней Олег Владимирович - ВГТУ, д-р физ.-мат.

наук, профессор, тел. +7(473)2466647

Костенко Виктор Викторович - ВГТУ, аспирант, тел. (473)

252-34-52

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

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

2. Математическая модель

2.1. Система уравнений

осесимметричного течения

Для определения параметров ламинарного неизотермического напорного течения

несжимаемой неньютоновской жидкости в цилиндрической трубе уравнения

гидродинамики и энергии для осесимметричного нестационарного потока записываются в безразмерном виде следующим образом: дп + 2 д(гу) дх г дг

0

(1)

du du „ du dp — + u — + 2v— = —— + dt dx dr dx

11 d , du

+ ~---d"(rB(2Peo ^ +

Re r dr dr

(2)

+

1

-----—)) + — — (2B-------------

Pe„ dx Re dx 2Pe

1 du

dx

)

0

ду ду

ду

---+ и-+ 2у— — —2Реп---+

дР

дt дх

1

дг

дг

1 д

дуч

+ — 2Ре0----------(г 2 В —) +

Яе

г дг

дг

(3)

1 д ди 1 ду

+ — Ре0 —(В(— +----------- —))

Яе 0 дх дг 2Ре02 дх”

д9 д9 „ де — + и — + 2у— — дt дх дг

1 э°_е Ре0 дх2

- +

(4)

1 д де 2

+ 4------------------(г —) + 4Вг0 ВИ

0

г дг дг

где

и

Ре

и— —; V — у— и и

0 .

г

г — —; я

х — ^; Р — Р — р0-

2Ре0 Я

t —

Ы

2Ре0 Я

ри2

Т — Т

е — ср

Т —Т

0 ср

— 2

Вг„

Ц 0и

МТ, — Тср)'

Яе

риЯ

„ иё Ре — — а

-0 J-cpJ Ц 0

Выше используются следующие

обозначения: и, V - компоненты безразмерной скорости; P - безразмерное давление; 0 -безразмерная температура; Т0 - температура жидкости на входе; Тср - температура окружающей среды; Яе - критерий Рейнольдса; Рг - критерий Прандтля; Ре - критерий Пекле; Вг - критерий Бринкмана; В1 - критерий Био; 1 -коэффициент теплопроводности; т -динамическая вязкость; р - плотность; a -коэффициент температуропроводности; В -безразмерная квазиньютоновская вязкость; Я -радиус трубы.

2.2. Расчетный алгоритм

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

Р+и — 2Ри + Р-и 2Ре2

2 Ах2

( ' Аг >

1 + —

+

Аг

Р

г

— 2 Р + Р

1,1 Гг, 1—1

V V 1,1 J

Рг+1,1 — ^—1,1_ + — —1

— (5)

2Ах

Аг

где ^ и О определяются из уравнений импульса и включают конвективные и диффузионные члены. После нахождения давления рассчитываются значения компонентов скорости на следующем временном шаге:

ип+1 — ип +А^ — дР);

уп+1 — уп + А(О — 2Ре2 ^.

Для решения уравнения энергии используется схема ВВЦП [1]:

е П+1 — е + Аt

Ре^Ах2 (е '*и — 2е - + е-и *

+ -

г Аг

‘, .7

2 (г1,1+1еУ+1 — (г1,1+1 + г>,;—1)е +

+ г1,1—1е у-О + 4ВГ0 В1, ;И5 —

(6)

2Ах

■(в „.,,—в ,ч, уА1 (в „1+,—о.,—,)

Предполагая, что изменение температуры мгновенно следует за изменением поля

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

К —1

В — к (Т )И/т , где т - показатель степенного реологического уравнения.

Диссипативная функция И равна:

И —

Г,

2Ре0Аг

\ ( Л2

VРе0ги ) I 4Ре0Ах

V

2Аг

4Ре02Ах

/

*

*

2

и

2

+

+

+

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

2.3. Граничные условия

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

Рис. 1. Классификация расчетных узлов

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

Узлы . О Это узлы, соседние с входными скоростями, причем снизу, сверху, справа лежат расчетные узлы, а слева - входные скорости, не меняющиеся по времени. Для этих узлов уравнение Пуассона выглядит так:

Р

г +1, у

Р,у 2Ре;

2Дх2

(- 2Р, у + Р

+

Дг2

(

Р

г, у+1

( А ^

Дг 1 + —

V

у-1

1 Дг

ЛЛ

г+1, у

м»-1, у

2ДxДt

20;

+ -

у/у 0 - + -

2Дх

- О Дг

(7)

г, у+1

и У-1

Нетрудно заметить, что в отличие от внутренних узлов расчетной области здесь в левой части уравнения отсутствует член

(Р,у - Р-1,у У(2&х), а в правой части вместо Е-1, у/ (2Дх) присутствует и- у /(2ДxДt)

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

Р - 2Р + Р

4+1,у ^4,у ^ -4-1,у

2Дх2

+ -

2Ре;

(

Дг

(

-Р

1 -

2Дг

Л

+

г, у

+ Р

1 -

2Дг

г

V г,у уу

20, , О

Е+1,у Ег-1,у

2Дх

+

(8)

+ -

,у

,у-1

гг ■ Дг

,у

Узлы О Это узлы, лежащие на выходной границе течения. При исследовании транспортных трубопроводов это сечение, которым мы искусственно «обрезаем» область расчета, считая, что течение установилось и

Эи Эу 0 у

— = — = и. Учитывая это, можно записать:

Эх Эх

2Ре;

( ( Р

г, у+1

20, у + 0 г

( Дг

1 +

г ,у

V ,у

, у +1 - О

Л

- 2Р + Р

г, у ^ Гг, у-1

С *

1 Дг

г

V г,у У/

(9)

", у

где

Дг

О. . =

г,у Яе

8рео г,у+1В,у+1У,у+1 - 2г,Д,уУ,.у + г,у-1В,.у-У,у-1

>7 г,у 1 ^ 1 ^-1 1 ^-1 г,у-

т, Дг2

,у

- Т7 (у-М+1 - у,у-1)

Узлы Н<. Эти расчетные узлы лежат на оси симметрии трубы и потока, которая, как

Эи 0

известно, характеризуется тем, что — = и и

Эг

у = 0, кроме того, конечноразностный аналог

г

Эу = Уу+1 - (~ уі,;+1) = уі, і+1 дг 2Аг Ат

замечаний можно записать:

. С учетом этих

Р - 2 Р + Р

4+1. і ^Гі, і^Гі-1, і

2 Ах2

Р+1,і Р-1, і

+

4Ре

2 Г Р.

і, і+1 - Р

V г. і у

20

(10)

2Ах

+

і, і+1

Аг

где

2(ві+ииі+и - 2Ві,іиі,і + ві-ииі-1,і)

-Т; ; = ------------------------------------------- ■

і,і

ЯеРеоАх

-(иі+1, і - “і-1,/ )

2Дх 4 г+1,у

Форма записи уравнений Пуассона для других узлов производится из подобных соображений. Необходимо отметить, что при дискретизации производных от давления таким образом сумма коэффициентов всей системы получающихся алгебраических уравнений равна нулю, что обеспечивает ограниченность расчетной схемы и выполнение интегрального уравнения баланса. Однако в некоторых расчетных узлах не обеспечивается диагональное преобладание матрицы коэффициентов, что обуславливает применение нейросетевого метода взвешенных невязок (НМВН) [4] для решения уравнения Пуассона.

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

к '(Тст - Тср ) =

ср' а Дг

где коэффициент теплопередачи от

внутренней поверхности трубы в окружающую среду с учетом толщины трубы и

теплоизоляционного покрытия равен

1

к'

1 + 5с а 2

+

Таким образом. перед

“ст '"из

каждым итерационным шагом ВВЦП определяется температура стенки по известной температуре в пристенке и заданному числу Б1:

А = 9ж

А /-«Т

1 + 1Б1Аг 2

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

соответственно, граничных условий. Граничные условия будут иметь особенности лишь на оси симметрии потока. на выходной границе и в точке пересечения этих двух классов. На оси симметрии схема ВВЦП упрощается (производные температуры по радиусу равны

нулю), либо приравнивается к нулю только —,

Эг

—Г г

дг V дг

Ч,і+1

Г е і і

Аг

і,і+1

Аг

На выходной границе можно ставить либо

— = 0, однако это условие будет формально Эх

означать, что профиль температуры установился и не будет меняться в дальнейшем. Это возможно лишь при установившемся профиле скоростей, при отсутствии вниз по потоку теплообмена со средой и если диссипация несущественна. Лучше

Э9

границе

выходной

разностей, а

Эх

представлять на с помощью левых

э 2е

Эх 2

0 , что формально означает

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

З Э9

границы. Здесь производная — представляется

Эх

левой разностью

э 2е

дг 2

как на выходной границе,

д 2 е

Эх 2

эе

0 и — = 0. дг

о

и

а

а

3. Нейросетевой алгоритм решения

уравнения пуассона

3.1. Постановка задачи обучения

нейросетевого пробного решения

Итак, на плоскости ХЯ определена цилиндрическая область, а в ней двумерная цилиндрическая сетка, заданная декартовым произведением двух одномерных сеток {хк},к = 1,...,п и {гI},I = 1,...,т . Под решением уравнения Пуассона будем понимать

нейросетевую функцию р = /нет х, г) , доставляющую минимум суммарной

квадратической невязки уравнения Пуассона, вычисленного по сеточным формулам типа (710) в совокупности узлов расчетной сетки. Представим пробное решение уравнения Пуассона в форме ч

p(w, х г) = £ у,/с (Ь + ^х + wi2г) + Ър . (11)

,=1

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

Для получения параметров решения уравнения (11), удовлетворяющих уравнению Пуассона, используются стандартные оптимизационные процедуры. Негативной особенностью НМВН является

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

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

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

Задача получения корректных параметров решения (11) наиболее важна на первом шаге реализации метода МАС. Здесь для задания начальных приближений весов и порогов сети могут использоваться результаты расчета распределения давления, выполненные с использованием классических вычислительных процедур. При расчете на следующих итерациях оптимизационный алгоритм настройки параметров нейросетевого решения (11) быстро находит решение, т.к. здесь по сути дела происходит подстройка распределения давления с учетом незначительного изменения правой части уравнения (5).

3.2. Сходимость вычислительного алгоритма

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

погрешностей, которые могут привести к рассогласованию и неустойчивости

приближенного решения, является дискретная аппроксимация производных, входящих в решаемые дифференциальные уравнения. Собственно, если использовать для нахождения параметров пробной функции аналитические решения интеграла взвешенных невязок, например, в одномерном случае | Вт (х)Я( х)ёх,

X

где Вт (х) - весовые функции, то в

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

множество точек, т.е. расчетная сетка конечных размеров. Очевидно, что если количество расчетных точек п ® ¥, то и интеграл

взвешенной невязки и сумма квадратов невязок

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

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

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

Условие устойчивости вычислительного алгоритма базового метода МАС

--------£ 0,25 в нашем случае было

Яе АхАг

модифицировано и на основании

многочисленных вычислительных

экспериментов принято следующее условие:

РеА1 £ 0,1.

Яе АхАг

3.3. Нейросетевой расчет распределения вихревой вязкости

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

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

Например, при повороте потока в канале на 90 градусов при Яе порядка 200 метод установления не дает сходящегося решения для гидродинамических параметров; в потоке возникают турбулентные флуктуации, приводящие к расходимости итерационного процесса. На рис. 2 приведено расчетное

распределение турбулентной вязкости в области течения согласно к-е модели, после определения которого используется тот же

гидродинамический алгоритм с заменой в

11

уравнениях движения на ——+ ц г.

---- на

Яе

решаются

Яе

совместно

обратного

Уравнения для кие нейросетевым методом

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

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

Рис. 2. Распределение турбулентной вязкости

4. Заключение

По известным распределениям скорости и температуры вычисляются параметры теплообмена и сопротивления в каждом сечении

трубы: градиент температуры и сопротивления в каждом сечении трубы, градиент температуры на стенке е'(0) , средняя температура жидкости е , местное число № = —2е'(0)/ е, среднее

— 2е'(0) „ 1 ф

число Ми =-------=---, местный С/ =---------— и

0

2Pe dx

среднии

1 X

Cf = — JCf (x)dx коэффициенты

x 0

сопротивления. Применение нейросетевых вычислительных технологий позволяет

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

неньютоновских степенных жидкостей.

Литература

1. Флетчер К. Вычислительные методы в

динамике жидкости: В 2-х томах.-М., «Мир», 1991.

2. Фройштетер Г. Б. Течение и теплообмен

неньютоновских жидкостей в трубах. Киев, «Наукова Думка», 1990. 216 с.

3. Шульман З. П. Конвективный

тепломассоперенос реологически сложных жидкостей. М., «Энергия», 1973. 351 с.

4. Кретинин А.В. Метод взвешенных

невязок на базе нейросетевых пробных функций для моделирования задач гидродинамики // Сиб. журн. вычисл. математики. РАН. Сиб. отд-ние. 2006. Т. 9. № 1. С. 23-35.

5. Кретинин А.В. Использование

динамических расчетных сеток в нейросетевом методе взвешенных невязок для моделирования гидродинамических задач / А.В. Кретинин, Ю.А. Булыгин, В.А. Волгин, В.Н. Апасов//

Нейрокомпьютеры: разработка, применение. 2007. № 9. С. 33-39

6. Кретинин А.В. Решение уравнений Навье-Стокса нейросетевым методом взвешенных невязок на базе нейросетевых пробных функций // Системы управления и информационные технологии. 2005. № 2(19). С. 17-19

7. Учет эффектов вязкоупругости в пристеночной области потоков нефтяных магистралей /

В.В. Фалеев, А.В. Кретинин, Ю.А. Булыгин, В.Д. Гриценко // Вынужденная конвекция однофазной

жидкости: Труды РНКТ-3. М.: Изд-во МЭИ, 2002. Т.2.

С. 285-288

8. Рубинский В.Р., Хрисанфов С.П.,

Климов В.Ю., Кретинин А.В. Математическое моделирование и экспериментальные исследования горения кислородно-метанового топлива при соосноструйной подаче в камеру сгорания ЖРД// Известия вузов. Авиационная техника. 2010. № 1. С 54-57

9. Нейросетевой алгоритм решения уравнений гидродинамики / Косинов О.А., Кретинин А.В., Стогней В.Г., Лукъяненко В.И., Вестник Воронежского государственного технического

университета. 2008. Т. 4. №4. С. 108-112.

10. A.V. Kretinin, Yu.A. Bulygin and M.I.

Kirpichev, “Method of Weighted Residuals on the Base of Neuronet’s Approximations for Computer Simulation of Hydrodynamics Problems”, Proceedings of IEEE 6th International Conference on Computational Cybernetics ICCC 2008. Stara Lesna, Slovakia. 2008. p. 237-240

11. Kretinin A.V. et al, “Neuronet's method of

weighted residuals for computer simulation of hydrodynamics problems”, Proceedings of 2010 International Joint Conference on Neural Networks, p. presented at WCCI 2010, Barcelona.

Воронежский государственный технический университет

THE MAC METHOD MODIFICATION USING NEURONET’S ALGORITHM OF THE POISSON’S EQUATION DECISION TO APPLY FOR HYDRODYNAMICS PROBLEMS

S.G. Valyukhov, A.V.Kretinin, O.V. Stognya, V.V. Kostenko

The MAC method is modified with neuronet’s numerical decision of the Poisson’s equation. The algorithm is applied for decision of Navier-Stokes equations with non-newton’s statement of viscous stress

Key words: hydrodynamics, not Newtonian liquids, neuronetwork modeling