CC BY
139
УДК 678.027
СИНТЕЗ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ТЕЧЕНИЯ ВЯЗКОЙ НЕНЬЮТОНОВСКОЙ ЖИДКОСТИ В КАНАЛЕ СО СЛОЖНОЙ ГЕОМЕТРИЕЙ И ИДЕНТИФИКАЦИЯ ЕЕ ПАРАМЕТРОВ С.В. Лавров
Статья содержит математическое описание течения вязкой неньютоновской жидкости в канале со сложной геометрией. Для описания характера течения были использованы уравнения Навье-Стокса, неразрывности и энтальпии. Выполнена идентификация параметров модели, которая позволяет с достаточной степенью точности определять изменение скорости и давления в процессе течения материала
Ключевые слова: моделирование, течение, неньютоновская жидкость
С целью моделирования течения вязкой неньютоновской жидкости в канале со сложной геометрией была выбрана отечественная система моделирования движения жидкости и газа FlowVision. Данный программный комплекс обладает широкими возможностями при задании условий трехмерного перемещения среды, а также позволяет различными способами визуализировать полученные результаты.
При формовании материала равномерность процесса зависит от сочетания различных характеристик, таких как давление, создаваемое шнеком, и сопротивление, обусловленное формующим каналом.
Для корректного проведения процесса моделирования был разработан алгоритм проведения математической обработки экспериментальных данных. Начальный алгоритм проведения процесса моделирования состоит из следующих основных этапов:
■ на основании теоретических и экспериментальных исследований выбор конструкций для проведения численного моделирования;
■ создание геометрии расчетной области во внешнем графическом редакторе и импортирование ее в оболочку FlowVision;
■ задание математических моделей для расчетной области;
■ задание граничных условий;
■ задание исходной расчетной сетки и критериев ее адаптации по решению и по граничным условиям;
■ задание параметров методов расчета;
■ проведение расчета;
■ просмотр результатов расчета в графической форме и сохранение полученных данных в файлах.
Поскольку сама программа FlowVision не содержит средства задания геометрии расчетной области, то для графического построения формующего узла была использована программа Компас 3Б V12. Построенный в ней формующий узел для большей
Лавров Сергей Вячеславович - ВГТА, канд. техн. наук, ст. преподаватель, тел. 8-920-442-25-99, e-mail:
наглядности на рис. 1, а представлен с вырезом четверти. В дальнейшем сохраненная в формате *.stl геометрия была уже импортирована в оболочку Flow Vision (рис. 1, б).
а б
Рис. 1. Геометрия формующего узла: а - созданная в Компас 3D V12; б - импортированная во FlowVision
Отличия между созданной во внешней программе и импортированной геометриями формующего узла обусловлены требованиями самой программы FlowVision [1]. Они способствуют наиболее корректной расстановке граничных условий на входе и выходе материала из формующего канала. Поскольку деталь является осесимметричной, расчетная область задачи представляет собой половину детали, что позволяет значительно снизить время расчета.
Характеристики формующего узла: длина формующего канала 1,1 м; угол конусной части 30°, длина конусной части 0,5 м, длина мундштука 0,3 м, сечение на выходе - квадрат 0,08x0,08 м.
При решении разнообразных инженерных задач, связанных с формованием пищевых масс, наиболее эффективными оказываются методы реологии, которые базируются на применении единых уравнений, справедливых для всех процессов, происходящих в рабочих каналах оборудования.
Основное свойство течения вязких неньютоновских жидкостей определяется развитием пластической деформации, представляющей вязкое течение, связанное с необратимым перемещением молекул и их групп на расстояние, превышающее размеры самой молекулы. Скорость развития пластической деформации зависит от температуры. В условиях установившегося течения вязкие материалы обладают свойствами так называемых аномально-вязких, или неньютоновских жидкостей. Это означает, что при весьма малых напряжениях
сдвига реологические свойства вязких материалов характеризуются постоянной ньютоновской вязкостью [2, 3]. Поэтому на следующем этапе выбирали модель течения и решаемые уравнения в расчетной подобласти [1, 4].
Была выбрана модель, характеризующая движение жидкости при малых скоростях и числах Рейнольдса менее 2300. Далее эту модель называем «Ламинарная жидкость» (уравнения Навье-Стокса, неразрывности, энтальпии):
- уравнение Навье-Стокса д¥ Ур 1
--+ V(V 0 V) = —±- +—V(pVV),
dt Р Р
- уравнение неразрывности
VV = 0,
- уравнение энтальпии
— + V(Vh) = — V dt р
—Vh
с
\
(1)
(2)
(3)
где р - плотность, кг/м ; V - скорость, м/с; р - давление, Па; ^ - динамическая вязкость, Па-с; к - энтальпия; Дж/кг; с - удельная теплоемкость, Дж/(кг-К); X - коэффициент теплопроводности, Вт/(м-К).
В качестве опорных величин были приняты температура Т = 273 К и давление р = 101000 Па.
При моделировании течения вязкой жидкости в модели «Ламинарная жидкость» уравнение Навье-Стокса дополняли граничными условиями. Для облегчения процесса расстановки граничных условий и исключения постановки несовместимых граничных условий, они объединены в определенные группы - типы границ, соответствующие конкретному физическому процессу, происходящему на границе. Для модели «Ламинарная жидкость» были выбраны следующие граничные условия:
I - граничное условие на непроницаемой твердой поверхности, называемое условием «стенка», характеризуемое «прилипанием» материала к стенке канала:
Щг=0; (4)
II - граничное условие, называемое «вход/выход», определяющее условия входа вязкой жидкости в расчетную область и обуславливающее направление вектора скорости материала. При этом на входе задавали дополнительно подтип граничного условия, называемый «нормальная скорость», в котором направление принимали строго перпендикулярно поверхности задания граничного условия, а значение вектора задавали предварительно;
III - граничное условие «свободный выход», определяющее выход материала в область с атмосферным давлением:
Р\г=0. (5)
Кроме того, использовалось также комбинирование граничных условий для одних и тех же расчетных поверхностей, так на границе «вход/выход» задавались граничные условия, как для скорости, так и для давления, что обусловлено, в первую очередь, рассчитываемыми величинами и точностью
проведения моделирования. Для ускорения расчета необходимо ввести дополнительное граничное условие симметрии, означающее, что проектируемая деталь является симметричной относительно вертикальной плоскости, проходящей через ее центр: V\„=0, (6)
V
dn
= 0.
(7)
В программе FlowVision используется декар-товая система координат, поэтому для расчета задавались изменяемой объемной сеткой, состоящей из набора ячеек заданного размера. Для уменьшения времени генерации и адаптации сетки задавались ячейками относительно большого размера, а в наиболее важных областях производили объемную адаптацию первого и второго уровня. В выбранных областях, где необходимо провести наиболее точный расчет, исходная ячейка разделялась на несколько одинаковых ячеек меньшего размера. Затем каждая из этих ячеек также делилась на ячейки еще меньшего размера. Дробление ячеек происходило до заданного уровня. С целью повышения точности расчета применялась «скошенная» расчетная схема, используемая при моделировании закрученных течений и позволяющая получить наиболее точные численные значения исследуемой величины [1, 5].
Для численного решения базовых уравнений в FlowVision использовали конечно-объемный метод, получивший широкое применение при моделировании перемещения жидкостей и газов. Этот метод положительно зарекомендовал себя для проектирования перемещения жидкостей и газов в широком диапазоне изменения их теплофизических свойств. Он основан на консервативных схемах расчета нестационарных уравнений в частных производных, которые по сравнению с неконсервативными схемами дают решения, точно удовлетворяющие законам сохранения (в частности, уравнению неразрывности) [6, 7, 8]. Метод конечных объемов сводится к аппроксимации сплошной среды с бесконечным числом степеней свободы совокупностью элементов, имеющих конечное число степеней свободы. Затем между этими элементами каким-либо способом устанавливается взаимосвязь. Популярность метода объясняется простотой его физической интерпретации и математической формы. Использование ЭВМ позволяет получать решения многих сложных технических задач включая моделирование течения вязких неньютоновских жидкостей, а использование в программе FlowVision древовидной структуры адаптивной сетки, при которой каждая ячейка связана с сеткой нулевого уровня, позволяет сократить время на генерацию расчетной сетки и снизить требования к оперативной памяти ЭВМ. Это положительно сказывается на скорости проведения расчета. Кроме того, при необходимости, для решения возникающей системы линейных алгебраических уравнений может использоваться как неявный (более надежный), так и явный (быстрее работающий, но расходящийся при больших шагах по времени) вариант итерационного процесса. Метод
базируется на эйлеровом подходе к описанию движения жидкости [1], суть которого состоит в том, что различные скалярные и векторные величины рассматриваются как функции переменных Эйлера от времени и координат точки в неподвижной системе координат.
В FlowVision численное интегрирование уравнений по пространственным координатам проводится с использованием прямоугольной адаптивной локально измельченной сетки. Такой подход обеспечивает, с одной стороны, использование простой равномерной неадаптивной сетки при решении задач с относительно несложной геометрией. С другой стороны, появляется возможность при решении задач со сложной геометрией проводить адаптацию (подстройку) сетки к особенностям геометрии вблизи границ, а при решении задач с разрывными течениями адаптацию по значениям искомых функций и их градиентов.
Процедура локального измельчения в области адаптации предусматривает возможность последовательного деления, начиная с исходной, каждой предыдущей ячейки на восемь более мелких ячеек до обеспечения выполнения условия адаптации (т. е. для достижения заданной точности вычисления искомой величины). На рис. 2 представлена схема локального измельчения адаптивной сетки, которая использовалась при моделировании течения.
При адаптации локальной сетки был применен метод подсеточного разрешения геометрии, предназначенный для аппроксимации криволинейных границ на прямолинейной сетке. При том ячейки, через которые проходит криволинейная граница, расщеплялись на несколько ячеек, теряющих первоначальную форму параллелепипеда и превращающихся в многогранники произвольной формы. Вместе с тем, никаких упрощений в математической модели для таких ячеек нет, что позволяет получать достаточно точные результаты расчета.
Свободная грань ячейки
Твердая грань ячейки
Рис. 2. Схема локального измельчения адаптивной сетки
Предварительно задавали сетку уровня 0 (самую грубую сетку, представленную на рис. 3, а). Для этого во вкладке «начальная сетка» в дереве варианта задавали направления сетки по осям х, у и г равные 10:10:20 соответственно. Неравномерное распределение плоскостей по различным осям связано с направлением движения моделируемых жидкостей.
а б
Рис. 3. Метод подсеточного разрешения геометрии: а - криволинейная поверхность проходит через ячейки; б - расщепление ячеек границей
Относительно большее количество расчетных ячеек по оси х обусловлено необходимостью точного определения оптимальной длины формующего канала.
Для повышения точности моделирования процесса произвели измельчение сетки первого порядка вблизи тех поверхностей, которым соответствует граничное условие «стенка» (рис. 4). Это позволило оптимально привязать расчетную сетку к заданной геометрии формующего узла.
Рис. 4. Измельчение сетки первого порядка
Следующим шагом являлось задание параметров метода численного моделирования. В программе FlowVision используется метод конечных объемов для численного решения уравнений конвективно-диффузионного переноса, записываемых в общем виде как:
+ Ч(У/) = Ч( ПУ/) + 0,
(8)
где / -рассчитываемая переменная, У - скорость, П - коэффициент диффузии, 0 - источнико-вый член.
При методе конечных объемов уравнения типа (8) интегририровали по объему каждой г-й ячейки расчетной сетки и по отрезку времени:
(9)
I/йу , п+1 -$/йУ г П + ^=
У I = I У I = I т8г
= Ц ПЧ/сксИ + Ц 0с1усИ
ту
где Уг - объем ячейки, Si - поверхность ячейки, Г, Г+1 - моменты времени начала и конца шага по времени.
1П+Х = Г + т .
(10)
При адаптации расчетной сетки и подсеточ-ном разрешении геометрии ячейка имеет форму произвольного многогранника. Те из граней, которыми ячейка граничит с другими ячейками, называется «свободными». Площадь /-ой свободной грани
в г-й ячейке обозначим через . «Твердыми гранями» являются те грани ячейки, которые образованы границей, пересекающей ячейку. Площадь /-ой твердой грани в г-й ячейке обозначим через gJ/ . Основываясь на этих обозначениях получим выражение (9) в разностной форме:
V (ГГ1 - /П) + Е / + Е ои/ + 0 = 0, (11)
/ /
где 0 - объемный источник переменной; / Г - ос-редненное значение переменной по объему ячейки в момент времени Г:
г
(12)
Осредненные плотности потоков рассчитываемой переменной и О/ через соответствующие грани з/ и g■ за шаг по времени определяли по
выражениям:
^ =\(/V+ич/)Л ,,
г ^
О = 1 (/V + П(У/ )^ )Л
ё,
(13)
(14)
где индексом V обозначены значения величин на границе расчетной области, которая соответствует грани g■.
Комплекс ПЧ/ , соответствующий диффузионному потоку величины /, аппроксимировали в FlowVision вторым порядком точности по пространственной переменной.
Одной из наиболее сложных задач при компьютерном проектировании процесса является аппроксимация конвективного потока в рассчитываемой геометрии. В зависимости от необходимой точности расчета используется несколько схем восстановления (восстановление первого порядка, гладкое восстановление высокого порядка и ступенчатое восстановление высокого порядка). Все эти схемы основаны на восстановлении рассчитываемой переменной / из ее средних значений внутри ячейки расчетной сетки и переноса восстановленной жидкости по линиям тока жидкости. Кроме того, дополнительно используется так называемая «скошенная» схема, предназначенная для восстановления величин в ячейках при диагональном перемещении заданной среды. Именно эта схема и была выбрана для дальнейших исследований процесса перемещения материала в формующем узле. Это обусловлено тем, что в проектируемой области присутствуют участки со «скошенным» направлением движения материала относительно ячеек сетки. На рис. 5 пока-
/ (х) =
зано расположение средних значений исследуемых величин в ячейках и схема аппроксимации уравнения конвективного переноса «скошенной» схемой. Иначе эта схема может быть обозначена, как схема аппроксимации «с дополнительной точкой». Ее записываем в следующем виде:
[ /П + /хП (Х+1/2 - ХХ ПРи Х ^ Х 1/й + Лй (Х - Х,-1/2Х ПРи Х < Х+1/2 - 1 где I - расстояние от стороны Х+1/2 до дополнительной точки, /■ и /и - величины скорости на левой и правой границах конечного объема. Их значения определяются из разложения функции в ряд Тейлора:
/(х) = /(Х+1/2) + /Х Ах +
+1 /Х Ах2 +1 /ХХХХ Ах3 + 0(Ат 4), (16)
где
2 ‘
Ах = х - х,
3!'
і+1/2 •
(17)
Коэффициенты выражения (16) находятся с использованием определений осредненного значения функции і-го конечного объема, а также соседних ему і-1 и і+1 объемов.
В результате получается:
/п = 0,5 (/,+1 + /г) - ОД67 (/,+1 - 2/, + 2/,-1);
/і = 0,5 (/, + /,-1) + °Л67 (/,+1 -2/, + 2/,-1).
(18)
и • • /+1 9
/ '
Х,-і:
Хі+і:
Хі-і:
Х1+1 2
а б
Рис. 5. Схемы расположения: а - расположение средних величин в ячейках; б - схема аппроксимации уравнения конвективного переноса «скошенной» схемой
Для достижения монотонности реконструкции величины / и / ограничиваются средними значениями соседних конечных объемов:
1 /г < /гг < /г+1, ПРи /г+1 - /г ^ 0;
I/г > /гг > /г+1, ПРи /г+1 - /г < 0;
[/г-1 ^ /1г < /г, ПРи /г - /г-1 ^ 0;
1 /г-1 > /1г > /г, ПРи /г - /г-1 < 0
Расстояние I и производные /хг, и /Хц на левой и правой сторонах У1 равны:
і =
/хі,
/хг,
/г, - /, і
і -1
(19)
Для численного решения уравнения Навье-Стокса используем неявный алгоритм расщепления по физическим переменным. Для этого запишем уравнения для движущегося объема ^:
і
| рСУ = | рСУ,
у, П
I рУСУ - | рУСУ = -11 РСШ + П.
У Пг тБ
(20)
где Б - поверхность объема ^, У - поле скорости рассматриваемой жидкости, р - плотность, П - члены уравнения Навье-Стокса, описывающие силу тяжести, вязкостные напряжения и т. п.
Разностный аналог уравнения Навье-Стокса выглядит как:
уИ+1урИ+1 - | рпупйу = т р&п+\з + Б, (У). (21)
В этом уравнении неизвестными величинами являются Уп+1 и Рп+1. Добавим и вычтем в нем дополнительные члены следующим образом:
. (22)
(Уп+' + У - У)У рп+! - 1рпУпСУ =
п
= -т(Е Р5п+гз - е Рьпь+Е РьпЪ)+вг (У)
Б Ъ Ъ
Это уравнение расщепляется на два:
УУрп+А - | рпУпСУ = -еРЪпь + пг (У), (23)
п, ъ
(Угп+1 - У)Угрп+1 = -ЕРБп+'•? + ЕР£Ъ . (24)
Б Ъ
В (23) используется поле давления, взятое на предыдущем шаге по времени. Это векторное уравнение представляет собой три уравнения конвективно-диффузионного переноса для трех компонент скорости жидкости. Для дальнейшего решения оно расщепляется следующим образом:
у, =■
1
(
Ур
п+1
I УпСУ -тЕ РьпЪ
\
У = У+У П (У).
(25)
(26)
Для определения поля давления рассмотрим условие несжимаемой жидкости, из которого видно, что:
Е уп+'з = 0,
(27)
Б
г п+1
ТТ'п+1
где УБ - значение скорости на границах конечного объема У .
Чтобы получить выражения для УБп+' запишем аналог уравнения (24), полученного интегрированием уравнений Навье-Стокса по движущейся
грани объема ^. Для грани этого объема, которая совпадает с гранью Ъ при =4 и с з при 1=1п+1 выра-
жение для УБп+1 будет иметь вид:
(ЧРп+1)
р
п+1
(ЧРп+1)
Р
п +1
. (28)
Подставляя полученное УБ (27) получим:
п+1
в выражение
е—т^1)
Б Ре
З = ЕУбЗ + Е^п-г(ЧРп+1)
Б ' ^ рп+1 Б Ъ Ръ
Ъ . (29)
Ъ
п+1
Б Б ь 1Ъ
После нахождения поля давления Рп+1 из (24) вычисляется поле скорости УБп+'. Для предотвращения осцилляции поля давления на неразнесенных сетках во FlowVision в уравнении (28) вводится разность между представлением градиента давления вторым и четвертым порядком точности [1].
После задания параметров метода численного моделирования производили запуск программы на генерацию обновленной сетки и расчет варианта по заданным условиям. Сам расчет происходит без вмешательства пользователя, поэтому отсутствует необходимость детального описания этого процесса.
Во время расчета представляется возможным отражение получаемых значений в графическом виде (заливка, вектора, изолинии и т. д.), а также получение характеристик исследуемой величины и сохранение данных в файл с возможностью обработки результатов во внешней статистической или математической программы.
Для отражения информации, получаемой при моделировании процесса, использовалась графическая визуализация результатов: изображение распределения переменной в плоскости методом цветовой (тоновой) заливки. Результаты расчетов сохранялись в виде отдельных графических файлов.
Для визуального отображения изменения скорости и давления использовали заливку цветом в зависимости от ее расчетного значения в каждой точке.
В каждом отдельном графическом файле определенный цвет соответствует конкретному значению исследуемой величины.
Для получения данных о скорости и давлении использовали параметр «характеристики» и задавали сохранение численных данных во внешний текстовый файл. По полученным результатам строили двумерные графики изменения величины скорости и давления по длине канала. На кривых V = /(I) и р = /(I) отражены участки частичной стабилизации давления в зоне прессования I, полной стабилизации в зоне формования II и падения давления в зоне III.
Изменение скорости течения материала в формующем канале показано на рис. 6.
г
т
ь
Б
и
Ъ
у
Как видно, скорость течения материала достигает максимального значения в зоне калибровки, поскольку давление в этой зоне снижается (рис. 6. и 7).
Рис. 6. Изменение скорости течения материала в формующем канале
Изменение давления в формующем канале
показано на рис. 7.
Р,
Па
101000
X
^1 ' 1 I * ^*
179000
171200
1634-00
155600
147300
140000
132200
124400
УІ 116600 108300 101000
канале
Рис. 7. Изменение давления в формующем канале
Результаты верификации модели с экспериментальными данными (рис. 8) свидетельствуют о её адекватности.
Рис. 8. Верификация модели с экспериментальными данными
Литература
1. Система моделирования движения жидкости и газа Ио%гУІ8Іоп 2.3. Примеры решения типовых задач [Электронный] / «Тесис», 2006. - 174 с.
2. Исаев, А.И. Инженерный метод расчета течения полимеров в каналах некруглого сечения [Текст] / А. И. Исаев, К. Д. Вачагин, А. М. Набережнов. // ИФЖ, 1974. -Т. XXVII, № 2. - С. 310 - 316.
3. Тадмор, З. Теоретические основы переработки полимеров [Текст] / З. Тадмор, К. Гогос. - М.: Химия, 1984. - 632 с.
4. Кондратин, Т. В. Применение пакетов прикладных программ при изучении курсов механики жидкости и газа [Текст] / Т. В. Кондратин, Б. К. Ткаченко, М. В. Бе-резникова, А. П. Зуев. - М.: МФТИ, 2005. - 104 с.
5. Примеры решения тестовых задач в программном комплексе FlowVision версия 2003 [Электронный] / ООО «Тесис», 2003
6. Андерсон, Д. Вычислительная гидродинамика и теплообмен. Ч. 1. [Текст] / Д. Андерсон и др. - М.: Мир, 1990. - 382с.
7. Госмен, А.Д. Численные методы исследования течения вязкой жидкости [Текст] / А. Д. Госмен, В. М. Пан, А. К. Рингел и др. - М.: Мир, 1972. - 324 с.
8. Сигерлинд, Л. Применение метода конечных элементов [Текст] / Л. Сигерлинд. - М.: Мир, 1979. -392 с.
Воронежская государственная технологическая академия
SYNTHESIS THE MATHEMATICAL MODEL FLOW A VISCOUS NON-NEWTONIAN FLUID IN A CHANNEL WITH COMPLEX GE OMETRY AND IDENTIFICATION OF ITS PARAMETERS S.V. Lavrov
Article contains a mathematical description of viscous non-Newtonian fluid in a channel with complex geometry. For the description of character of a current the equations of Navier-Stokes, continuity, and enthalpy. Completed identification of model parameters, which allows a sufficient degree of accuracy to determine the change in velocity and pressure in the flow of material
Key words: modeling, flow, non-newtonian fluid
