CC BY
96
УДК 536.491:53.044:519.633.2
МЕТОД ФАЗОВОГО ПОЛЯ В ЗАДАЧЕ ОПИСАНИЯ ПРОСТРАНСТВЕННОЙ ДИНАМИКИ ПОВЕРХНОСТИ ПОЛИСТИРОЛА ПРИ ТЕРМИЧЕСКОЙ ДЕСТРУКЦИИ
1 ЛЕБЕДЕВ В.Г., 1 ЛЕБЕДЕВА А.А., 2ОБУХОВ А.А.
1Удмуртский государственный университет, 426034, г. Ижевск, ул. Университетская, 1 2Физико-технический институт Уральского отделения РАН, 426000, г. Ижевск, ул. Кирова, 132
АННОТАЦИЯ. В работе предложена фазово-полевая модель, описывающая пространственную динамику поверхности при деструкции полистирола под действием температурного фронта в бескислородной среде. Параметры фазово-полевой модели сопоставлены с параметрами кинетики процесса деструкции, экспериментально полученными на основе PCFC (pyrolysis combustion flow calorimetry) метода. Проведено моделирование одномерной фазово-полевой задачи, иллюстрирующее особенности численного расчета в рамках представленной модели.
КЛЮЧЕВЫЕ СЛОВА: динамика поверхности, пиролиз полистирола, фазовое поле. ВВЕДЕНИЕ
Описание динамики поверхностей в трехмерном пространстве является весьма нетривиальной задачей, если предполагать наличие резкой, постоянно меняющейся границы, разделяющей пространственные области с кардинально отличающимися свойствами [1]. Альтернатива состоит в использовании диффузной границы, предполагающей плавный переход от одной области к другой. Один из таких подходов связан с методом фазового поля [2]. Важное преимущество метода фазового поля, по сравнению с другими подходами к задаче описания движущихся поверхностей, состоит в возможности его термодинамического обоснования. Поэтому метод фазового поля давно и успешно используется при описании фазовых переходах первого рода [3], задачах образования трещин [4] и т.д. В частности, метод фазового поля применялся для описания перехода из метастабильного состояния в энергетически более выгодную фазу за счет химической активации [5]. Еще одним преимуществом метода фазового поля является простота обобщения алгоритма на произвольную размерность пространства. При этом усложняются формулы, и увеличивается объем вычислений, но сам принцип остается тем же самым.
ЛИТЬЕ ПО ВЫЖИГАЕМЫМ МОДЕЛЯМ
Типичная задача описания пространственных поверхностей возникает при компьютерном моделировании литья по выжигаемым моделям (Lost Foam технология -современный вид точного металлургического литья [6], занимающего прочные позиции на мировом рынке. Физический процесс по Lost Foam технологии сводится к изготовлению прототипа из вспененного полистирола (далее EPS - expanded polysterene) с последующим помещением прототипа в оболочку (формовочную смесь) и заливкой расплавленным металлом [7]. При соприкосновении с металлом EPS газифицируется и частично уходит в форму, частично поглощается металлом. Наличие оболочки ограничивает доступ кислорода, поэтому процесс деструкции полистирола не приводит к сгоранию продуктов распада.
Поглощение продуктов деструкции расплавом создает предпосылки для возникновения внутренней микропористости металла, что отрицательно сказывается на механических свойствах изделий, особенно из алюминия. Положительные стороны технологии - точная форма отливки и относительная дешевизна процесса. Использование EPS низкой плотности уменьшает количество газа и, тем самым, существенно снижает
проблемы с пористостью. Тем не менее, разработка математических моделей, описывающих процесс взаимодействия расплавленного металла с EPS, остается актуальным направлением в создании программных средств технологического прогнозирования литейных процессов [8].
Эксперименты с EPS показывают, что в силу быстроты процесса заливки, деструкция полимера происходит лишь в непосредственной близости от металла [9], рис. 1. Таким образом, в первом приближении можно выделить область расплавленного металла, газовую прослойку и EPS в твердой фазе. Между данными областями происходит интенсивный тепломассообмен, исследование которого требует знания пространственного расположения фронта деградации EPS пены.
В целях простоты, в данной работе мы ограничимся выводом фазово-полевой модели для пространственного описания поверхности выжигаемого полистирола и численным решением модельной задачи, иллюстрирующей метод фазового поля.
Рис. 1. Частичное заполнение формы металлом [9]. Слева затвердевший металл, справа - остаток EPS между слоями краски
КИНЕТИКА ПИРОЛИЗА ПОЛИСТИРОЛА В ОТСУТСТВИЕ КИСЛОРОДА
Тепловая деградация полистирола в бескислородной атмосфере, например, в азоте, изучалась во многих работах [10, 11] с помощью термогравиметрии и дифференциальной сканирующей калориметрии при скоростях нагрева 0,5^20 К/мин. С помощью безмодельного изоконверсионного метода была исследована зависимость энергии активации процесса деструкции от параметров. В результате исследований было обнаружено, что энергия активации является достаточно постоянной и равной примерно 200 кДж/моль в широком диапазоне конверсии. Несмотря на то, что многоступенчатый механизм часто встречается в полимерах, из полученных данных был сделан вывод, что деструкция полистирола происходит в одну стадию, приводя к появлению сухого остатка и горючих газов. Работа [12] на основе PCFC метода (pyrolysis combustion flow calorimetry) подтвердила ранее сделанные выводы, показав, что пиролиз полистирола в бескислородной атмосфере действительно можно формально рассматривать как одношаговый процесс. Поэтому в дальнейшем, мы не будем рассматривать процессы плавления полистирола, считая, что при достижении определенной температуры полистирол распадается из твердой фазы в горючую газовую смесь. Из рассматриваемых типов реакции в работе [12], наилучшее соответствие получено для автокаталитической реакции вида:
= A exp(-E / RT )(1 - c)09 (0,05 + Z06 ), (1)
где предэкспонента равна A = 1,18-Ю11 c-1, а энергия активации вычислена равной
m — m
E = 168 кДж/моль. Величина % =—0- характеризует оставшуюся массовую долю
m - m¥ 0¥
полистирола в момент времени при этом то - начальная масса полистирола, т - текущая масса полистирола в момент времени т- масса сухого остатка после деструкции. Вид кинетической кривой приведен на рис. 2.
Рис.
0,0 0,2 0,4 0,6 0,8 1,0
Массовая доля х, безразм.
2. Вид автокаталитической модельной функции [12]
МОДЕЛЬ ФАЗОВОГО ПОЛЯ ДЛЯ ДЕСТРУКЦИИ ПОЛИСТИРОЛА
Идея фазового поля основана на введении локального параметра порядка - фазового поля ф( г ,0, характеризующего состояние среды в точке пространства г в момент времени (. Для удобства обычно принимается, что ф меняется в интервале от 0 до 1, причем границам интервала сопоставляются различные состояния среды (например, разные фазы). Введение параметра порядка дает возможность термодинамического описания релаксационного изменения состояния среды [2] на основе уравнений типа Гинзбурга-Ландау [13]:
д(р = -м ар{(р,Т ]
(2)
дГ д(р
где М- кинетический коэффициент (мобильность), Е[ф,Т] - свободная энергия системы, управляющая процессом релаксации.
В случае деструкции полистирола под двумя состояниями будем понимать исходное состояние недеградировавшего полистирола и конечное состояние, когда от полистирола остается некоторый сухой остаток и горючие летучие фракции. Относительно состава газов можно полагать, что он является равновесным, поскольку скорость химреакций существенно превышает характерные скорости процесса выжигания полистирола расплавленным металлом. Но в данной работе вопрос о составе и кинетике формирования химического состава продуктов деструкции рассматриваться не будет.
Для процесса распада полистирола естественно интерпретировать фазовое поле как массовую долю, заменив для удобства х = 1-ф. Тогда исходному состоянию полистирола соответствует ф = 1, а конечному ф = 0. Соответственно, кинетическое уравнение (1) перейдет в уравнение
( = А ехр(-Е / ЯТ) ((0,05 + (1 - ()0,6). (3)
д(
Для описания термодинамического состояния деструктирующего полистирола в качестве свободной энергии выберем интеграл по области занятой полистиролом:
( = Д 2 е2У(2 + Wg (()
йУ,
где градиентное слагаемое учитывает наличие границы [14], коэффициент г определен ниже, g(ф) - модельная функция [2] вида:
g () = 4(2(1 -()2.
(4) будет
(5)
Функция g(<p), заданная "double-well" потенциалом (5), соответствует наличию двух равновесных состояний полистирола с потенциальным барьером W между ними. В равновесном случае SF/бф = 0 получаем уравнение:
где А - оператор Лапласа,
Aj - Wg'(j) = 0, (6)
dg 1
gj) = j(1 -j)(1 - 2j). dj 2
Для одномерной задачи с граничными условиями
lim j(z) = 0 и limj(z) = 1,
это уравнение имеет хорошо известное решение в виде "кинка":
1L J z - z w
2
d
jo(z) = - 1 + th—± I, (7)
описывающее равновесие двух состояний среды по направлению z (z0 - граница) c полушириной фронта 6, где:
d = V2e2 / W. (8)
Вклад градиентного слагаемого в выражении (4), вычисленный с помощью решения ф0, обычно отождествляется с поверхностной энергией [2]. В случае полистирола, поверхностную энергию E0 можно оценить как энергию активации слоя толщины 26:
г pEz+d pEd
E0 =- \j0(x)dx °-,
m z»-d m
2 2 3
где p - плотность полистирола р ~ 1,5-10- ^ 6,5-10- г/см ), / - молярная масса, среднее значение которой для полистирола можно принять /л ~ 950,4 кг/моль [12]. Приравнивая
1 +¥
E0 = {jx))2dx ° E0, окончательно определяем параметры модели
W = 12pE e = pr (9)
m m
Чтобы граница пришла в движение, необходимо нарушить равновесие, сделав состояние, соответствующее исходному полистиролу, метастабильным, рис. 3. В рамках фазового поля, метастабильное состояние можно получить с помощью непрерывного преобразования плотности свободной энергии F[$,T], добавив вклад движущей силы (driving force) перехода fFD, умноженный на "размазанную" по толщине диффузной границы функцию "ступеньки". Нетрудно видеть, что в нашем случае роль движущей силы играет выражение
fFD(T)=F0exp(-E/^T), (10)
а в качестве "ступеньки" можно взять интеграл от оставшейся части кинетической функции:
р( j) = N\js0,9 (0,05 + (1 - s)0,6 )ds, (11)
J0
где N = 3,51 - нормировочный коэффициент, вычисленный из условия p(1) = 1.
График функции р(ф) представлен на рис. 4. График свободной энергии с метастабильным состоянием на рис. 3 получен линейной комбинацией функций g(p) и р(ф).
Рис. 3. Изменение свободной энергии за счет движущей силы приводит к появлению метастабильного состояния
Рис. 4. Вид функции "ступеньки" р(ф)
Таким образом, для описания метастабильного состояния при деструкции полистирола, к свободной энергии (4) достаточно добавить интеграл по объему:
Ffd [j, T ] = J fFD (T) p(j)dV.
(12)
Тогда, из уравнений Гинзбурга-Ландау (2) следует уравнение фазового поля, описывающее динамику границ при деструкции полистирола:
j = M[в2Dj- Wg'(j)-F0 exp(-E/RT)p'(j)) (13)
Физический смысл уравнения достаточно прозрачен. При RT<<E источник fDF экспоненциально стремится к нулю, поэтому в системе устанавливается равновесие с выбранной шириной границы. При повышении температуры, величина силы fDF(T)p}(y) растет, в результате чего "кинк" приходит в движение. Форма "кинка" определяется конкуренцией диффузии (Дф) и вклада "double-well" потенциала (g'(^)).
Величина мобильности M обычно выбирается на основании экспериментов по движению границы. Однако, в данной задаче такое сопоставление с экспериментом возможно лишь после полного описания процесса, учитывающего, в том числе, движение свободной границы расплавленного металла, выход и распространение продуктов деструкции. Не имея такой возможности сопоставления, отметим, что в представленной модели фазового поля отсутствует собственный характерный масштаб времени. Действительно, скорость деструкции существенно зависит от температуры нагрева и функция fFD(T) является быстро растущей функцией вплоть до температур T~104. Поэтому в качестве характерного масштаба времени возьмем время t0 теплового процесса и перейдем к безразмерному времени t ^ t0t. Рассматривая фазовое поле как средство для отслеживания пространственных характеристик, выберем для простоты
m
toM
6pE
(14)
Учитывая, что для установившегося профиля фронта
£2А(-Wg() » 0,
движущая сила должна воспроизводить кинетическую функцию (2), находим
A
t oMFo
10 A / N.
(15)
что дает окончательное уравнение для динамики границ при деструкции полистирола в виде:
( = 52А( - 2g' (() - А0 ехр(- Е / ЯТ)р ((). (16)
W
ВЫДЕЛЕНИЕ ТЕПЛА НА ГРАНИЦЕ ДЕСТРУКЦИИ
Процессы перераспределения тепла вблизи диффузной границы выжигаемого полистирола подчиняются закону сохранения внутренней энергии:
— + У • У = 0, (17)
Эt v ;
где и - плотность внутренней энергии, У = —1Т - поток тепла, X - коэффициент теплопроводности. В силу термодинамического соотношения между плотностями равновесных термодинамических потенциалов
и = I — т 1 ,
J эт'
производная по времени от плотности внутренней энергии может быть записана через плотность свободной энергии 1(Т,ф) как
Эи =_т Э21 т+[ I—т Л
эт эт2 ^ Эр ЭтЭр
р.
Пользуясь определением, для единицы объема имеем теплоемкость:
С =—т
С т Эт2'
с помощью которой, из закона сохранения (17), получаем уравнение теплопроводности в виде:
Су Эт = У(IVт)_ Г1 — т ^Цр. (18)
Эt
Величина
Эр ЭтЭр
б
г Э1 Э2/ ^
— т
Эр ЭтЭр
характеризует мощность тепловыделения в единице объема за счет перехода из метастабильного в устойчивое термодинамическое состояние. В качестве равновесного вклада в свободную энергию системы при таком феноменологическом подходе необходимо учесть плотность свободной энергии вида
1о = ^(р) + /ро (т)р'(р).
Подставляя эту плотность в источник, с учетом (14) и (15), находим
б = ^о ехр(—Е / Ят) • р '(р)(Е / Ят — 1) — Wg '(р). (19)
Выбирая характерный пространственный размер Ь, а также характерные значения теплопроводности С0 и теплопроводности Х0, перейдем к безразмерным переменным г ^ гЬ, X ^ Х0Х, Су ^ С0СУ, что определяет характерное время распространения тепла ^ = С0Ь}/Х0 и мощность тепловыделения, пересчитанную в градусы Кельвина:
б = Г (А ехр(—Е / Ят) • р '(р)(Е / Ят — 1) — 2g '(р)). (20)
Здесь слагаемое с функцией g'(ф) описывает изменение внутренней энергии при переходе через потенциальный барьер, поэтому интеграл по времени от данного вклада в фиксированной точке пространства будет равен нулю:
| g Хр)рЛ = g(р)|10 = 0
хотя это слагаемое влияет на динамику процесса во времени. Слагаемое с функцией р'(ф) учитывает изменение внутренней энергии между метастабильным и устойчивым состоянием. Поскольку в интервале температур Т < 2000 К множитель (Е/ЯТ-1) всегда отрицателен, р'(ф) > 0, и ( < 0, вклад' в тепловыделение от движущей силы~р'(ф) всегда положителен, в отличие от знакопеременного вклада g'(ф).
МОДЕЛЬНЫЙ РАСЧЕТ
В качестве иллюстрирующего примера, рассмотрим решение одномерной задачи для уравнений (16) и (18). За единицу длины выберем величину S. Используя имеющиеся визуальные данные по деструкции, находим из рис. 1, что характерная величина S может быть оценена как величина порядка S ~ 1 мм.
Теплофизические характеристики будем считать постоянными в каждой среде. В расчетах использованы следующие теплофизические данные для EPS: теплопроводность ^-EPS = 0,035 Вт/(м-К), теплоемкость CEPS = 1,45 кДж/(кг-К), плотность pEPS = 0,02 г/см , взятые из [8]. В качестве теплофизических характеристик газовой фазы использованы данные по метану (при 1000 К) [15]: теплопроводность ANH4 = 0,037 Вт/(м-К), теплоемкость CNH4 = 4,475 кДж/(кг-К), плотность pNH4 = 1,96^ 10-4 г/см3. Характерное время прогрева в газовой фазе на расстоянии S = 1 мм получается равным t0 = CNH45 /ANH4 ~ 0,025 c. Вычисление константы A0 в соотношении (15) дает значение A0 ~ 8,4-10 .
Пусть пространственная область соответствует интервалу 0 < z < 100^. Будем полагать, что область 0 < z < 10^ не меняется со временем и заполнена расплавленным металлом с температурой T0 = 1000 К. Выбор T0 соответствует расплаву чистого алюминия, перегретому относительно точки плавления Tm = 933,3 К. Область 10^ < z <z1 - газовая фаза, которая, в данной модельной задаче, обеспечивает перенос тепла за счет эффективной теплопроводности и никак не влияет на процесс деструкции. Область z1 < z < 100^ заполнена твердым полистиролом. Пространственное разбиение выбрано равным N = 2^ 10 . Масштаб по времени определяется из условия, чтобы температуропроводность в газовой фазе равнялась единице. Тогда безразмерная температуропроводность полистирола равна 0,032.
Граничные условия для температуры заданы в виде условий I рода T(0) = T0, T(100^) = Tp = 300 К. Начальное распределение температуры выбрано в виде:
T0, z < 10d, T( z,0) = (Tp, z >10d.
Для моделирования уравнения теплопроводности (18) используется абсолютно неявная разностная схема, решаемая методом прогонки [16]. Для моделирования уравнения фазового поля (16) применяется простейшая явная схема по времени. Устойчивость схемы проверялась эмпирически. Расчеты показывают, что устойчивость выполняется при Дt < 0,45Дх . Выбор
явной схемы обусловлен двумя обстоятельствами. Во-первых, в силу сильной нелинейности уравнения фазового поля, расчет с большим шагом по времени приводит к большим погрешностям в решении. Поэтому, введение неявных алгоритмов с целью увеличения шага по времени и сокращению времени расчета, в такой задаче теряет смысл. Во-вторых, использование неявных алгоритмов при решении уравнений фазового поля приводит к рассогласованию между диффузионным вкладом и "double-well" потенциалом, в результате чего меняется профиль фазового поля и, соответственно, транспортные характеристики границы между различными состояниями.
В качестве граничных условий для фазового поля использовались условия вида:
= 0,
z=100d
а в качестве начального условия - решение (7) при z0 = 10^.
j(0) = 1; j
dz
Отметим также одну из особенностей численного решения уравнения для фазового поля в рамках предложенной модели. Она связана с наличием ненулевых производных у используемой кинетической функции на концах интервала, рис. 2. Строго говоря, это означает, что точки ф = 0 и ф = 1 не являются неподвижными точками для уравнения фазового поля при численном решении задачи. Поэтому, при моделировании вместо
функции g(ф) (5) использовался потенциал со стенками
~<„)=^0 <Р< 1
[ ¥, ф< 0 и Р > 1.
Результаты моделирования представлены на рис. 5. Смещение фронта деструкции показано на рис. 6. Из представленных графиков видно, что по мере распространения тепла фронт деструкции продвигается совместно с тепловой волной, сохраняя свой профиль.
Смещение фронта соответствует корневой зависимости от времени ~ -Л . При приближении тепловой волны к границе 2 = 100^, фронт фазового поля несущественно замедляется. Из графика на рис. 5 также видно, что тепловыделение на фронте (20) не вносит существенного вклада в профиль температуры на фронте деструкции полистирола.
Рис. 5. Профили температуры и фазового поля в различные моменты времени
Рис. 6. Смещение фронта деструктировавшего полистирола от начального положения с течением времени
ЗАКЛЮЧЕНИЕ
В работе показана возможность построения фазово-полевой модели, динамика которой имитирует процесс деструкции полистирола при литье по выжигаемым моделям (Lost Foam технология). В основу динамики фазового поля заложены данные по пиролизу полистирола, экспериментально полученные PCFC-методом. Движение границы деструкции полистирола в рамках предложенной модели качественно проиллюстрировано на примере численного решения одномерной задачи и легко обобщается на случай трехмерного пространства.
Для полного описания процессов взаимодействия расплавленного металла с полистиролом при литье по выжигаемым моделям представленная фазово-полевая модель деструкции полистирола должна быть дополнена моделями, описывающими движение поверхности расплава и взаимодействие горячих продуктов разложения полистирола с расплавом и формой. Разработка соответствующих моделей является необходимым шагом для создания программных продуктов, моделирующих технологические процессы литья по выжигаемым моделям.
СПИСОК ЛИТЕРАТУРЫ
1. Gurtin M.E., Voorhees P.W. The thermodynamics of evolving interfaces far from equilibrium // Acta Mater. 1996. V. 44, № 1. P. 235-247.
2. Provatas N., Elder К. Phase-Field Methods in Materials Science and Engineering. Weinheim : Wiley-VCH, 2010. 254 p.
3. Galenko P., Jou D. Diffuse-interface model for rapid phase transformations in nonequilibrium systems // Phys. Rev. 2005. E 71. P. 046125.
4. ^rma A., Levine H., ^ssler D. Phase-Field Model of Mode-III Dynamic Fracture // Phys. Rev. Lett. 2001. V. 8704. P. 5501-5504.
5. Provatas N., Grant M., Elder К. Phase-field model for activated reaction fronts // Phys. Rev. 1996. E 53. Р. 6263- 6272.
6. URL: http://www.achfoam.com/Lost-Foam-Casting-Process.aspx (дата обращения 20.11.2014).
7. Григорьев В.М. Литье по выжигаемым моделям. Хабаровск : Изд-во ХГТУ, 2002. 57 c.
8. URL: http://wp lvm.mkmsoft.ru/ (дата обращения 20.11.2014).
9. Sen I. Degradation Mechanism of Expanded Polystyrene (EPS) Foam in Lost Foam Casting, PIPS Approach for Synthesis and Novel Expansion Techniques for Cellular Foam. PhD diss., University of Tennessee, 2007. 145 p.
URL: http://etd.utk.edu/2007/SenIndraneel.pdf (дата обращения 5.12.2014).
10. Peterson J.D., Vyazovkin S., Wight C.A. ^netics of the thermal and thermooxidative degradation of polystyrene, polyethylene and poly(propylene) // Macromol. Chem. Phys. 2001. V. 202. P. 775-784.
11. Chen К., Harris К., Vyazovkin S. Tacticity as a factor contributing to the thermal stability of polystyrene // Macromol. Chem. Phys. 2007. V. 208. P. 2525-2532.
12. Snegirev A.Yu., Talalov V.A., Stepanov V.V., Harris J.N. Formal kinetics of polystyrene pyrolysis in non-oxidizing atmosphere // Thermochim. Acta. 2012. V. 548. P. 17-26.
13. Гинзбург В.И., Ландау Л. Д. К теории сверхпроводимости // Журнал экспериментальной и теоретической физики. 1950. Т. 20. C. 1064-1081.
14. Cahn J.W., Hilliard J.E. Free energy of a nonuniform system. I. Interfacial free energy // J. Chem. Phys. 1958. V. 28. P. 258-267.
15. Таблицы физических величин : Справочник / под ред. акад. И.К. Кикоина. М. : Атомиздат, 1976. 1008 с.
16. Самарский А. А. Теория разностных схем. М. : Наука, 1977. 656 с.
PHASE-FIELD DESCRIPTION OF THE SPATIAL DYNAMICS OF A POLYSTYRENE SURFACE DURING HEAT DESTRUCTION
1Lebedev V.G.,1Lebedeva A.A., 2Obuhov A.A.
:Udmurt State University, Izhevsk, Russia 2
Physical-Technical Institute, Ural Branch of the Russian Academy of Sciences, Izhevsk, Russia
SUMMARY. A phase-field model is proposed for description of spatial dynamics of polystyrene surface under temperature front influence in oxygen-free environment. Parameters of phase-field model are derived from the parameters of the degradation kinetics, which was experimentally obtained on the basis of PCFC (pyrolysis combustion flow calorimetry) method. Simulation of one-dimensional phase-field problem is done. Calculations illustrate the features of numerical modeling in the framework of this model.
KEYWORDS: dynamics of surface, pyrolysis of polystyrene, phase field.
Лебедев Владимир Геннадьевич, кандидат физико-математических наук,
заведующий кафедрой теоретической физики УдГУ, тел. (3412)916-130, e-mail: lvg@,udsu.ru
Лебедева Алла Анатольевна, методист УМР физико-энергетического факультета УдГУ, e-mail: allaleb1985@gmail. com
Обухов Александр Андреевич, аспирант ФТИ УрО РАН, e-mail: Obuh87@bk.ru
