CC BY
13
УДК 519.6:532.62
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПОЛЗУЩЕГО ТЕЧЕНИЯ РАСПЛАВОВ ПОЛИМЕРНЫХ СИСТЕМ. Часть 2. Алгоритмы численного решения
АЛЬЕС М.Ю.
Институт механики Уральского отделения РАН, 426067, г. Ижевск, ул. Т. Барамзиной, 34
АННОТАЦИЯ. Рассматриваются эволюционные схемы численного решения краевых задач ползущего неньютоновского течения реакционно-способных расплавов полимерных систем.
КЛЮЧЕВЫЕ СЛОВА: полимерные неньютоновские расплавы, реокинетика, ползущее течение, краевые задачи, метод конечных элементов, эволюционные численные схемы.
РАСЧЕТ ПОЛЕЙ ТЕМПЕРАТУРЫ И СТЕПЕНИ ПОЛИМЕРИЗАЦИИ
Для области V , ограниченной поверхностью £ = 8Т и , запишем слабую форму Галеркина для уравнения энергии (см. (15) в [1]) с естественными граничными условиями УгТ| = Я + а(Т - Тср) на Бч
я
\lVTVyJV + \[рср (Э? + ЦУгТ) - а]^ + \[я + а(Т - Тср )\¥я<&, (1)
V V V
и уравнения макрокинетики реакции полимеризации (см. (16) в [1])
- Е_
{[э = | к/^уЬУпГР, (2)
V V
при некотором известном распределении гидродинамических параметров ц, I . ,т.
(ц - скорость, I. - тензор-градиенты скоростей, т. - напряжения сдвига). Здесь
£Т, - части границы, на которых задаются соответственно температура и закон
теплообмена с окружающей средой; я - плотность теплового потока; а - коэффициент
теплообмена; Тср - температура окружающей среды; Э^ - частная производная по (•).
Заметим, что, так как диффузией степени полимеризации пренебрегается (см. (16) в [1]), то перенос ¡3 осуществляется только за счёт движения жидкости. Поэтому самостоятельных
граничных условий для 3 не выставляется. В начальный момент времени t0 задаётся
распределение температуры Т0, степени отверждения ¡0, скорости Ц, положение
свободной поверхности (исходная область У0). В данной постановке предполагается, что
превалирующим фактором, влияющим на процесс течения, является изменение эффективной вязкости за счет: 1) реализуемых тензор-градиентов скорости течения; 2) изменения температуры и степени полимеризации. Объемные изменения в движущейся среде в связи с изменением температуры и степени полимеризации в данной постановке не учитываются.
Произведя триангуляцию р и представив элементные распределения температуры и
степени полимеризации в виде ( Гп - множество узлов элемента п )
Т = УТ , Р = ¥Д , £ГЯ, (3)
будем иметь две взаимосвязанные системы сеточных дифференциальных уравнений относительно неизвестных узловых значений температуры
^ ат | рсръту1 +2 ^ | т [1(У у хуу+
п=1 V п=1 V
у п у п
N
+ Рсрь' (Уу ] dVn + 2 лп аД | [аТу + ч - аТсР ] уд<£ -
N
-2 ат/[еуэд+диуг^)+дт]у^п = о, (4)
п=1 V
и степени полимеризации
N N -
2ат\[у ЭЛ +ргиугуг:\¥^п = 2а т|кре *Ту(Р)у^п (5)
п=1 V п=1
у(Д) = (1 -уЛ)(1 -cyW.fi,), (т = 1,M;1 ,у,кеГпДеВп) ,
где N, М - число конечных элементов и узлов в триангуляции р; Гп - множество узлов элемента; X - множество элементов, у которых хотя бы одна грань используется при аппроксимации £; Вп - множество узлов элемента, находящихся на £; а 7т - логическая процедура
а г
т
1, если узел У е Г конечного элемента п
Г1, "п е Ь
соответствует узлу т = 1,М триангуляции р ЛпЬ = <
I 0, в противоположном случае.
0 в противоположном случае;
Символ е обозначает, что при суммировании по впереди стоящему немому индексу производится последовательный перебор элементов данного множества; индексом п отличается принадлежность п-му конечному элементу. Граничные условия на £ являются
естественными для уравнения энергии. Граничные условия на £Т : Т|5 = Т* удовлетворяются
путём непосредственной замены уравнения для соответствующего узла, например т, тождеством Тт = Т*. Здесь £Т - части поверхности £, на которых задается распределение
Т*
температуры Тт .
Запишем уравнения (4),(5) в виде
сТ+к1Т+^ - ^(Д Л) = о, (6)
с2ЛЛ = 02Д, Т), (7)
где точка обозначает производную по времени. Воспользовавшись неявной схемой для (6),(7) будем иметь две взаимосвязанные системы алгебраических уравнений относительно температуры Т+А( и степени полимеризации Д+м на временном шаге ^ + &
(с = сТ +М[01(Д^,Д -^], (8)
сД^ = сД + Ь02Д+*, Т^). (9)
В целях сокращения записи, временная идентификационная метка I + & у операторов К, G1, ^ опущена.
Сеточная задача (8), (9) является нелинейной. Простейшим приемом, при помощи которого может быть решена проблема нелинейности является вычисление правой части в (9) по решению, полученному в предыдущий момент времени
= сЬ + Лгв2(&, Т ). (10)
Однако отметим, что такой приём эквивалентен применению явной схемы для аппроксимации уравнения макрокинетики (7), которая является лишь условно устойчивой и накладывает очень жесткие ограничения на шаг Лt [2, 3]. Воспользуемся методом простой итерации
х+1 ^+1
(с +ЛК1) Т^ = сТ +Л^(Д+^, Д) - Р ], (11)
х+1 х х
С2 Д+Л = сЛ +ЛtG2(bt+лt,). (12)
Из условия локальной сходимости метода для уравнения макрокинетики, обуславливающего нелинейность задачи, имеем следующую оценку допустимого шага по времени в системе (11), (12) (в "вмороженных" координатах)
< ехр(^Г) . (13)
К & [1 + Со(1 - 2Д+Л)]
Здесь отметим, что в рассматриваемых задачах шаг Лt во многих случаях целесообразно определять исходя из кинематического условия на свободной поверхности и выставлять в (8), (9) в качестве заданного параметра, определяемого гидродинамикой процесса, геометрией стенок. В этом случае удовлетворение условию (13), как показали расчёты, как правило, оказывается обременительным с точки зрения вычислительных затрат. Невыполнение же условия приводит к возникновению неустойчивости сначала для полей /3 , а затем расходится всё решение. Для преодоления указанных трудностей введём вместо сеточного аналога уравнения макрокинетики на слое t + Лt (9) следующее эволюционное уравнение
&+Л! + сА* = сЛ +ЛtG2(Л+лt, Т+Л ), (14)
х = 0 =л,
где X/ - фиктивное безразмерное время.
Аппроксимируя (14) по явной схеме с шагом ЛХ/ > 0 , запишем следующий итерационный алгоритм решения задачи (8), (9)
х+1 х+1
(с +ЛК1) Т^ = сТ +&[01(Я+Л< Л) - Р ], (15)
х+1 х х х х
Ь t+* = Ь t+Лt -ЛХ [с (А+* -Л) - лЮ2 (Д+*, Т<+лt)]. (16)
Из условия локальной сходимости разностного аналога эволюционного уравнения для Ьг+Лг (в "вмороженных" координатах) имеем следующую оценку области допустимых значений итерационного параметра ЛХ/ , обеспечивающего сходимость процесса (15), (16)
2
ЛХ <--• (17)
1 + ЛК& ехр - —— [1 + со(1 -2Д+*)]
V у
Численные расчёты показали, что алгоритм (15), (16) при выполнении условия (17), которое удовлетворяется без каких-либо практических трудностей, обеспечивает всегда устойчивое сходящееся решение.
Рассмотрим теперь вопрос об аппроксимации конвективных членов. При использовании симметричных функций для аппроксимации конвективных членов в уравнении энергии и макрокинетики получаемое численное решение, как показали расчёты, часто является неустойчивым в областях, прилегающих к границам области, на которые натекает среда. Для получения устойчивого решения известно несколько способов. Самым простым (в некоторых случаях единственным) является измельчение расчётной сетки конечных элементов. Другим способом является применение схем против потока. Применительно к методу конечных элементов схема против потока может быть реализована путём представления конвективного члена производной по направлению и использования аппроксимаций типа конечно-разностных [4, 5]. В качестве направления выбирается направление, противоположное вектору скорости в данной расчётной точке области. Такой метод, однако, приводит к появлению дополнительных узлов и сложным аппроксимационным выражениям.
Другой способ реализации схемы против потока в методе конечных элементов заключается в использовании тензора коэффициента диффузии (вязкости, теплопроводности), который имеет несколько большие значения в направлении против вектора скорости [6]. Такой способ лишён эффекта появления большой схемной вязкости, присущего схемам против потока, но он достаточно сложен в реализации.
При использовании метода взвешенных невязок могут быть также использованы весовые функции, имеющие смещение в направлении против потока [4,6]. Такой способ введения схемы против потока, как и перечисленные выше, приводит к несимметричной матрице коэффициентов конечно элементных уравнений.
Воспользуемся следующей аппроксимацией [7]. При вычислении вкладов конвективных членов в систему проекционно-сеточных уравнений для узла конечно-элементной сетки вместо значений Т ( , Д принимаются значения Ти,Ди в точке, расположенной на расстоянии в направлении, противоположном узловой скорости
и . При этом, в качестве Т|,Д| принимаются значения, полученные на предыдущем временном слое. Такой способ в совокупности с измельчением сетки конечных элементов, как показали расчёты, позволяет получать устойчивое решение. Расчёты показали также, что в областях значительных градиентов полей Т и Д схемы против потока сглаживают получаемое решение. В этих случаях получение удовлетворительного результата возможно только путём измельчения сетки.
РЕШЕНИЕ СВЯЗАННЫХ ЗАДАЧ ГИДРОРЕОДИНАМИКИ И ТЕПЛООБМЕНА
Решение задачи гидродинамики реакционного формования полимерных изделий сводится к решению двух взаимосвязанных систем нелинейных сеточных уравнений относительно скорости иг и давления р (см. (10),(11) в [8]), температуры Т и степени полимеризации Д (4),(5). Распределение температурно-конверсионных полей зависит от реализуемой картины течения. С другой стороны, температура и степень полимеризации существенно влияют на эффективную вязкость полимерной массы.
Решение гидродинамической части задачи при фиксированных распределениях температуры и степени полимеризации обсуждалось в [8]. Выше были рассмотрены способы получения устойчивых сходящихся решений уравнений энергии и макрокинетики реакций полимеризации при неизменном распределении скорости по области интегрирования.
Рассмотрим возможные алгоритмы совместного решения нелинейной системы гидродинамических уравнений и нелинейной системы температурно-конверсионных уравнений в момент времени I + Лt . Совместное решение систем существенно усложняется тем, что положение свободной поверхности, а, следовательно, и область интегрирования, в момент времени t + Лt является так же неизвестным задачи. Определение положения свободной поверхности на данном временном шаге с помощью какого-либо итерационного процесса в совокупности с итерационными схемами определения других неизвестных приводит к чрезвычайно сложным и трудоёмким алгоритмам решения задачи. Реализация таких алгоритмов сопряжена с очень большими затратами счётного времени. Практически реальным, с этой точки зрения, является определение формы свободной поверхности по известному распределению скорости с предыдущего шага по времени. Тогда решение систем проекционно-сеточных уравнений будет осуществляться при известном положении свободной поверхности. Форма свободной границы аппроксимируется сторонами прилегающих к ней конечных элементов, её новое положение будет определяться следующим кинематическим соотношением
X1 (г + Лt) = хи ^) +1 Ц ^ , I = 1, £, (18)
где £ - число узлов, лежащих на свободной границе.
Системы гидродинамических и температурно-конверсионных уравнений являются нелинейными и для решения каждой из них используются эволюционные численные схемы, изложенные в [8] и выше (см. (14)). При совместном решении этих систем возможно объединение итераций по нелинейностям каждой системы в один итерационный процесс следующим образом. При известном начальном приближении температуры, степени полимеризации и скорости производится одна итерация для уравнений энергии и макрокинетики. По полученным распределениям температуры и степени полимеризации вычисляются реологические параметры жидкости и отыскивается приближённое распределение скорости и давления на первой итерации по нелинейности уравнения движения и удовлетворения условию несжимаемости. Затем производится следующая итерация в системе температурно-конверсионных уравнений и т.д.
Однако, нелинейность гидродинамических уравнений не связана явно с решением уравнений энергии и макрокинетики, а нелинейность последних также не связана явно с решением задачи гидродинамики. Поэтому такой алгоритм, как показали расчёты, имеет медленную сходимость при неудачном начальном приближении.
Более эффективным является следующий алгоритм. Начальное приближение температуры и степени полимеризации используется для определения значений реологических параметров и проводится цикл итераций для уравнения движения и условия несжимаемости. После удовлетворения условиям точности по скорости и давлению выполняется цикл итераций по нелинейности температурно-конверсионной системы уравнений с использованием полученного приближённого распределения скорости. Полученное поле температуры и степени полимеризации используется для поправки реологических параметров и выполняются итерационные циклы следующего приближения.
При достаточно малых шагах Лt, как показали расчёты, можно без потери в точности выполнять только по одному циклу итераций описанного выше алгоритма на данном временном слое, если для решения одной из проекционно-сеточных систем использовать решение другой системы с предыдущего шага по времени. Как правило, расчётный шаг по времени определяется гидродинамикой процесса и геометрией стенок из соотношения (18) при условии допустимого перемещения фронта свободной поверхности на величину, не превышающую минимальный по области шаг конечно-элементной сетки. Получаемые при этом значения шагов по времени являются много меньшими характерных времен тепловых тс и химических тр взаимодействий Лt <<тс ~ ]ЗрсрЯ~1 , Лt <<тр ~ ехр(Е / ЯТ)(КдС0)-1.
Кроме того, на каждом шаге по времени решается стационарное уравнение движения. Поэтому без потери точности, как показали расчёты, можно построить следующую схему вычисления скорости, давления, температуры, степени полимеризации на данном временном шаге.
Сначала решается система уравнений движения и неразрывности, реологические параметры в которой вычисляются по значениям температуры и степени полимеризации, снесенным с предыдущего временного шага. После достижения требуемой точности в итерационном процессе решения нелинейного уравнения движения и выполнения условия несжимаемости, полученное распределение скорости используется для решения систем нелинейных сеточных уравнений энергии и кинетики полимеризации. После этого решение задачи на данном временном шаге считается найденным.
СПИСОК ЛИТЕРАТУРЫ
1. Альес М.Ю. Математическое моделирование ползущего течения расплавов полимерных систем. Часть 1. постановка задачи // Химическая физика и мезоскопия. 2015. Т. 17, № 2. С. 208-213.
2. Рихтмайер Р., Мортон К. Разностные методы решения краевых задач. М. : Мир, 1972. 420 с.
3. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М. : Энергоатомиздат, 1984. 152 с.
4. Thomasset F. Implementation of Finite Element Methods for Navier-Stokes Equation. New York : Springer-Veklag, 1981. 159 р.
5. Bristeau M.O., Glowinsky R., Periaux J., Perrier P., Pironneau O. A Finite Element Approximation of Navier-Stokes Equations for Incompressible Viscous Fluids // Proceedings of IUTAM Symposium: Springer-Verlag. 1980. Р. 78.
6. Leonard B.P. A Stable and Accurate Modelling Procedure Based on Quadratic Upstream Interpolation // Computer Methods in Applied Mechanics and Engineering. 1979. V. 19(1). Р. 59-98.
7. Bonque J.P., Ibler B., Labadie G. A Finite Element Method for Navier-Stokes Equations // Proceedings of Int. Conf. Numer. Meth. in Non-Linear Problems. 1980. V. 1. P. 709.
8. Альес М.Ю. Эволюционные схемы численного решения нелинейных краевых задач неньютоновских ползущих течений // Химическая физика и мезоскопия. 2013. Т. 15, № 2. С. 193-200.
MATHEMATICAL MODELING CREEPING MELT FLOW POLYMER SYSTEMS. Part 2. Algorithms for numerical solution
Alies M.Yu.
Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, Izhevsk, Russia
SUMMARY. We consider the phenomenological descriptions and differential formulation of the problem creeping non-Newtonian flow of reactive melt polymer systems.
KEYWORDS: non-Newtonian polymer melts reokinetika, creeping flow, boundary value problems, finite element method, numerical evolutionary scheme.
Альес Михаил Юрьевич, доктор физико-математических наук, профессор, главный научный сотрудник ИМ УрО РАН, тел. 89128563824, e-mail: [email protected]
