Численное исследование динамики фазового разделения в бинарной полимерной смеси при спинодальном распаде на больших временах Текст научной статьи по специальности «Физика»

ВЫСОКОМОЛЕКУЛЯРНЫЕ СОЕДИНЕНИЯ, Серия А, 1997. том 39, № 6, с. 1014-1020

ФАЗОВЫЕ - ПРЕВРАЩЕНИЯ

УДК 541.64:536.7

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ДИНАМИКИ ФАЗОВОГО РАЗДЕЛЕНИЯ В БИНАРНОЙ ПОЛИМЕРНОЙ СМЕСИ ПРИ СПИНОДАЛЬНОМ РАСПАДЕ НА БОЛЬШИХ ВРЕМЕНАХ

© 1997 г. Е. В. Простомолотова*, И. Я. Ерухимович*, JI. И. Маневич**

* Московский государственный университет им. М.В. Ломоносова. Физический факультет

117234 Москва, Воробьевы горы ** Институт химической физики им. H.H. Семенова Российской академии наук

117977 Москва, ул. Косыгина, 4 Поступила в редакцию 11.07.96 г.

Принята в печать 18.09.96 г.

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

Исследование динамики фазового разделения бинарных полимерных смесей представляется важным как с теоретической, так и с практической точек зрения. Действительно, диффузия макромолекул характеризуется большими временами релаксации, что облегчает наблюдение образующихся структур [1-3]. Кроме того, для описания фазового разделения в полимерных смесях развита эффективная процедура в рамках теории среднего поля. Наконец, характерные свойства таких полимерных систем позволяют создавать новые уникальные материалы, находящие широкое применение на практике и открывающие новые горизонты в развитии современных технологий.

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

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

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

Описание процесса спинодального распада в линеаризованном случае

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

^ = (11У(А?га(1(ц/ЛГ)) (1)

ш

Здесь ф(г) - локальное значение объемной доли компоненты А, г — время, Т- температура, к - постоянная Больцмана, ц(г) = 5^&р(г) - локальный химический потенциал, равный вариационной производной свободной энергии системы F по

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ДИНАМИКИ ФАЗОВОГО РАЗДЕЛЕНИЯ

1015

<р(г) и А - коэффициент Онзагера, который определяется выражением [9]

Здесь

Л =

^аО-ФО^ + ^ВФО*8'

(1а)

^аР = -^^(П'СФо)*

+ 2

2,2 а к

36ф0(1

5ф(г, 0 = £с(к)е

-1/1(к) /кг

е ,

х-'(к) = * Л(ф0)^

,2\

1-

к2сг = -К

п18фо(1-фо)

(7)

где тА, тв - микроскопические времена релаксации, //А и //в - степень полимеризации компонент, Ые - число мономеров между зацеплениями вдоль цепи, а - длина статистического сегмента, предполагаемая для простоты одинаковой у обоих компонент.

В рамках решеточной модели Флори свободная энергия ^для бинарной полимерной смеси записывается в виде

Р_ = гГРр а\Уф)2 ^

кТ J Цг 36ф( 1 - ф)У ^

^ = $^ + (1-ф)1п(1-ф) + х(р(1_ф) (3)

(% - параметр взаимодействия Флори-Хаггинса). На ранних стадиях развития процесса, когда отклонение локального значения объемной доли ф(г) от ее среднего значения ф0 еще мало (6ф 1, 8ф = ф-фо), уравнение (1) можно линеаризовать относительно 5ф и решить его аналитически, что было впервые сделано Каном, Хиллардом [5] (см. также работы [6-9]).

Повторим вкратце их основные результаты. Пусть в пространстве создано некое начальное распределение концентраций ф(г, 0). Представим функцию 5ф(г, 0) = ф(г, 0) - фо в виде разложения

Фурье 5ф(г, 0) = УАс(к) еАг и подставим его в линеаризованное относительно 5ф уравнение (1). Тогда это уравнение в к-представлении примет вид

(4)

где ^о - вторая производная от Решением уравнения (4) являются волны

(5)

экспоненциально растущие или убывающие в зависимости от знака декремента тг'ОО. который определяется следующим выражением:

(6)

Из соотношения (6) следует, что при ГЦ /кТ > 0

(к2сг < 0) система находится выше спинодали, и все волны 5ф(к) затухают с течением времени. Если же ГЦ /кТ < 0 (кгсг > 0), то амплитуды волн с \к\ < ксг с течением времени экспоненциально нарастают. При этом, как следует из формулы (6), наиболее

быстро растут волны с к„ = к„142, которые по амплитуде обгоняют все остальные, и на начальном этапе процесса формируются домены с характерным масштабом к = 2п/кт, который называется длиной кановской волны.

Численное решение уравнения диффузии на поздних стадиях эволюции в одномерном и трехмерном случаях

На более поздних стадиях спинодал ьного распада отклонение ф от ф0 уже перестает быть малым, и линеаризация, проводимая на начальном этапе, оказывается неприменимой. В этом случае необходимо решать нелинейное уравнение диффузии, что как правило можно сделать лишь численными методами. Подобные расчеты проводились для решения уравнения (1) в одномерном [10, 11] И трехмерном [12, 13] случаях, которые показали возможность существования двух различных сценариев поведения полимерной системы на больших временах.

Первый сценарий, наиболее ярко наблюдавшийся при решении одномерного уравнения [11], характеризуется тем, что переход системы к состоянию полного фазового разделения осуществляется через образование промежуточных "кинетически стабильных" (т.е. долгоживущих) пространственно неоднородных структур. Между этими долгоживущими структурами происходят достаточно быстрые переходы, связанные с укрупнением масштаба структур, причем времена этих переходов на порядок меньше времени жизни кинетически стабильных структур. Для трехмерного случая такой сценарий имеет место лишь при небольших охлаждениях, а при глубоких охлаждениях реализуется второй сценарий, когда наблюдается скейлинговое поведение системы и характерный масштаб структуры плавно возрастает со временем по закону ~*а (а = 1/4—1/3 [12,13]).

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

ся интересным исследование поведения таких систем при решении уравнения (1) в двумерном случае.

Численное решение уравнения спинодального

распада на поздних стадиях эволюции в двумерном случае и проблема естественного начального распределения (р(г)

Уравнение (1) решалось на квадратной области со стороной квадрата, кратной длине канов-ской волны, для симметричной бинарной полимерной смеси с ф0 = 0.5 и = Рассматривалось два вида начального распределения: регулярное и случайное. Первое из распределений имеет вид

(Р(г,0) = Фо + е5т(^8т^, (8)

/

где п,т- целые числа; L -Lia- обезразмеренная

длина кановской волны, равная L = 132.46; р -число кановских волн, укладывающихся на стороне квадрата рассматриваемой области. Такое начальное условие позволяет рассмотреть эволюцию отдельной плоской волны, близкой к волне максимального роста.

В результате численного решения уравнения (1) (с начальным распределением вида (8), где мы положили р, п, m - 4) были получены профили функции 6ф(г, г), структурного фактора

5(k,i) = Jexp(ik(r,-r2))8<p(r„0x

х6ф(г2, t)dVidV2/V

(9)

L

0.010

0.006

0.002

20

60 t

100

и корреляционной функции G(r, t) =

= а также построены зависимости

свободной энергии F системы и скорости процесса v(t), которая определялась как норма в ¿2 производной Эф Idt выражением

t 2 1/2

v(i) = (J((ç(r,f + Ai)-9(r,f))/Af)2dr)

На первой стадии процесса в случае регулярного начального распределения на временах порядка времени роста кановской моды наблюдается увеличение амплитуды начального распределения, затем (при t - 6) рост замедляется и образуется кинетически стабильная долгоживу-щая структура, за время существования которой профиль концентрации практически не меняется. Это хорошо видно на зависимости свободной энергии от времени (рис. 1), где на том же промежутке времени свободная энергия также практически не меняется до t « 85, когда симметрия начального распределения нарушается и начинает резко меняться вплоть до нового квазисгабильно-

Рис. 1. Зависимости свободной энергии ДО системы (сплошная линия), где f[t) = 105F(f), и скорости процесса v(r) (штриховая линия) в случае регулярного распределения, (t - в единицах t(*m))-

го состояния. Хорошо иллюстрирует это и зависимость скорости процесса v(t) = Эф/Э/, на которой долгоживущим структурам отвечает нулевая скорость, а резким переходам - увеличение скорости (рис. 1).

Другой тип информации о структурной эволюции системы представлен на рис. 2, где, в частности, показано изменение структурного фактора S(k) системы при переходах от одной квазисгабиль-ной структуры к другой. Видно, что вместо четырех максимумов в структурном факторе S(k, t), относящихся к начальному синусоидальному распределению, появляются два, соответствующие ламелярной структуре.

Очевидно, однако, что регулярное начальное распределение вида (8) практически невозможно создать в реальной полимерной смеси. На практике начальный профиль плотности является случайным, причем характер распределения возможных начальных функций 0ф(г) можно определить из следующих соображений. Пусть система достаточно долго выдерживалась при некоторой температуре Т0 выше спинодали, после чего мы быстро (практически мгновенно в момент времени t = 0) охладили ее до температуры Г ниже спинодали. Очевидно, что при этом начальный профиль плотности 0ф(г) отличается от константы лишь в силу термических флуктуаций плотности, характер которых зависит от значения температуры Т0, при которой приготовляется полимерная система. Другими словами, поскольку скорость релаксации состава системы значительно меньше скорости релаксации ее температуры, в ней сохраняется первоначальное распределение плотностей как "замороженное воспоминание" о начальном состоянии при высоких температурах.

Рис. 2. Эволюция различных характеристик системы в ходе спинодального распада для регулярного начального распределения вида (8). Сверху вниз изображены численные значения функции ф(г), структурного фактора 5(к) и корреляционной функции С(г) для двумерных областей (х, у) и (кх, ку). Величина этих функций отображается оттенками серого цвета, г = 4 (а), 84 (б), 88 (в) и 99 (г).

Поэтому начальный профиль плотности естественно строить в соответствии с теорией флук-туаций [13], согласно которой амплитуды разложения 8ф(г) в ряд Фурье распределены по Гауссу с дисперсией

1Ы> ■ ¿У

1

МАФО МВ(1-ФО)

2.2 а к

-2Хо +

(Ю)

18ф0(1

:_1

-Фо)]

Здесь V- объем системы и 2%0 = 9/Г0, где 6 - температура Флори и Т0 - температура, при которой система была приведена в термодинамическое равновесие перед последующим охлаждением до температуры Т ниже спинодали. Для дискретной равномерной сетки начальное распределение имеет вид

(11)

Ф(Г) = Фо + Аф(г),

где

N/2

п,т = -N/1

^[Ъфгх + щ0 + 2яч); рЬ

соображениям. Поскольку реальные системы всегда трехмерны, в качестве начального распределения возьмем двумерное сечение трехмерного распределения в объеме V = (рЬ)2з (я - толщина слоя). При этом дисперсию (10) необходимо просуммировать по одной из координат кг. Заменим суммирование интегрированием. При этом считая, что на дискретной сетке, состоящей из NxNyNz узлов (Ых = УУу = АО,

а\к2х + к2у + к2) = <*ЕЩ£} + ЩИ, (12) {р1) "

получим формулу

Х<|Дф*|2> = -Цх

>2.2

(р1)3

N,12

¿1

(13)

1__2% + +А'2)

"дФо ЛГв(1-Фо) (/71)218фо(1_фо)

Т| — случайная величина, равномерно распределенная на сегменте [-0.5; 0.5]; Ы- число узлов сетки, на которой решалось уравнение; ^ = Ъ,(кх, ку) - случайная величина, распределенная по Гауссу с дисперсией, вычисляемой согласно следующим

N,/2

(Р^1/2 1

¿1

(2тс )У + ]2 + 12)

Фо ВД-Фо) (Д)218ф0(1-ф0)

(а)

(в)

ш ш Ш ш

■ ■ ■ ■ ■ j V 1_

• у г щг " щ ь â ф Ж * аЯ > # Я * m Ш 1 Ф JtL H ™ # ■ — 1 D v m » » m -W-.

Рис. 3. Эволюция различных характеристик системы в ходе спинодального распада для случайного начального распределения (11). г = 0 (а), 4 (б) и 12 (в).

В результате выражение примет вид £<|Дф*|2> = о(Хо)-о(0),

где функция о(Хо) определена следующим образом:

(14)

<*(Хо) =

л/П

;-2Хо +

2%{р1) ^аФо ^в(1-ф0) ВЫСОКОМОЛЕКУЛЯРНЫЕ СОЕДИНЕНИЯ Серия А том 39 № 6 1997

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ДИНАМИКИ ФАЗОВОГО РАЗДЕЛЕНИЯ

1019

„ 2,.2 .2. Л 8л: 0 +J )

^ 9 (Р1)

-1/2

агОД

(2тсЛГ/*)

18х

х

1

МдФо ^в(1"ФО)

-2Хо +

8 п(12 + /)" 9 (р1)

-1/2 \

Обратим внимание на то, что если к неограниче-но, что имеет место в непрерывном случае, когда /V, —► оо, I = х, у, г, ряд (11) расходится. Даже при конечных не очень больших значениях N1 (соответствующих в наших расчетах 16 узлам на одну кановскую волну), дисперсия <р(г) может оказаться порядка единицы, что приводит к появлению нефизических значений <р > 1 или <р < 0. Именно для преодоления подобных эффектов из дисперсии вида (14) с Хо * 0 и вычиталось выражение того же вида с % = 0 (см. формулы (13), (14)). (Это соответствует обрезанию высоких частот в начальном распределении.) Отметим также, что как уменьшение толщины слоя 5, так и ее увеличение (вплоть до бесконечности) при фиксированном сеточном шаге А = не влияет на значение интеграла (14) и определяется только зависимостью от А. Поэтому без ограничения общности

можно положить 51 = рЬ и Nz = N, т.е. рассматривать двумерное сечение трехмерного распределения в кубе. В случае же уменьшения сеточного шага А (вплоть до среднего расстояния между соседними мономерами а) дисперсия (14) будет расти, поскольку флуктуации в ячейке с характер-

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

Для качественного выяснения роли начального флуктуационного распределения состава мы решали уравнение (1) для шести различных случайных реализаций начального распределения, для каждого из которых вычислялись те же характеристики, что и в описанном выше случае регулярного начального профиля плотности (8). Для анализа роли размера системы рассматривались области двух размеров: со сторонами квадрата,

равными 2Ь и 4£, что соответствует N = 32 и 64 сеточным узлам. Исходные данные для расчета брались следующие: Ые = 50; = Л/в = 1000; % =

= 0.0025; Хо = 0.0005; фо = 0.5; тА = тв = 10-11 с; I =

= 132.46. Сеточный шаг Л = Ь /Ы.

Как и в предыдущем случае, на начальном этапе эволюции происходит быстрое убывание высоких гармоник с к > кт и рост гармоник с волновым числом к < кт, т.е. преобладают наиболее быстрорастущие (кановские) волны. На рис. 3 показано состояние шести систем для двух моментов времени / = 0 и 12. Видно, что с течением времени структуры укрупняются, амплитуды гармоник с волновыми числами к < кт нарастают, а высшие гармоники (с к > кт) убывают (во втором ряду рисунков от облака пиков, существующего при 1 = 0, остается всего 2-6 пиков).

В действительности образовавшаяся структура характеризуется масштабом порядка канов-ской волны уже на временах I ~ 1-4. Когда ампли-

0.016

0.012

0.008 -

0.004

Рис. 4. Зависимость средних квадратов амплитуд отдельных гармоник от времени. Кривые 7, /', 2, 2', 3, 3\ 4, 4' соответствует гармоникам с волновым вектором к = (0, 4), (4, -1), (1, -3), (3, -1), (2, -1), (2, -2), (0,2), (1, -2).

Рис. 5. Зависимости свободной энергии от времени для шести систем (в области 64 х 64 узлов) со случайным начальным распределением и, для сравнения, системы с регулярным начальным распределением (см. рис. 1).

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

что для области размером 4 L х 4 L с течением времени растут более низкие гармоники с волновыми числами (0,2) и (1,-2), а гармоники с канов-ским волновым числом (0,4) спадают. Видно, что с течением времени начинают преобладать все более низкие гармоники (например, ¿о2 вместо

/~2 2

&о4), где мы ввели обозначение k,j = km*Ji + j ¡p.

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

На рис. 5 изображены зависимости свободной энергии F(t) от времени для шести систем со случайным распределением вида (11) и, для сравнения, системы с регулярным распределением вида (8) для двумерной области со стороной квадрата,

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

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

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Hashimoto T., Такепака M., Izumitani Т. // J. Chem. Phys. 1992. V. 97. С. 679.

2. Такепака М., Izumitani Т., Hashimoto Т. // J. Chem. Phys. 1993. V. 98. С. 3528.

3. Lauger J., Lay R., Gronski W. // J. Chem. Phys. 1994. V. 101. C. 7181.

4. Лифшиц E.M., ПитаевскийЛ.П. Физическая кинетика. M.: Наука, 1979. Гл. 12.

5. Cahn J.W., Hilliard J.E. //J. Chem. Phys. 1958. V. 28. C. 258.

6. Cahn J W. Ц Acta Metall. 1962. V. 10. С. 179.

7. Cahn J.W. H Acta Metall. 1966. V. 14. С. 1685.

8. Cook Н Е. II Acta Metall. 1970. V. 18. С. 297.

9. De Gennes P.G. II J. Chem. Phys. 1980. V. 72. C. 4756.

10. Митлин B.C., Маневич Л.И., Ерухимович И.Я. H Журн. эксперим. и теорет. физики. 1985. Т. 8. № 12. С. 495.

11. Митлин B.C., Маневич Л.И. И Высокомолек. соед. А. 1988. Т. 30. № 1.С. 9.

12. Kotnis М.А., Muthukumar M. // Macromolecules. 1992. V. 25. С. 1716.

13. Brown G., Chacrabarti A. //J. Chem. Phys. 1993. V. 98. № 3. C. 2451.

14. Ландау Л.Д., Лифшиц E.M. Статистическая физика. M.: Наука, 1976. Ч. 1. Гл. 12.

A Numerical Study of the Long-Time Dynamics of Phase Separation during Spinodal Decay in a Binary Polymer Blend

E. V. Prostomolotova*, I. Ya. Erukhimovich*, and L. I. Manevich**

* Faculty of Physics, Lomonosov Moscow State University, Vorob'evy Gory, Moscow, 117234 Russia ** Semenov Institute of Chemical Physics, Russian Academy of Sciences, ul. Kosygina 4, Moscow, 117977 Russia

Abstract—The process of two-dimensional spinodal decay in a binary polymer blend was studied by numerically solving the corresponding nonlinear diffiision equation with a preset initial inhomogeneity. Evolution of the free energy, composition distribution, and structural factor was analyzed for both a regular initial distribution (with the Fourier expansion containing a small number of harmonics) and a random distribution (modeling typical inhomogeneities occurring in the system as a result of thermal fluctuations). It is demonstrated that the formation of kinetically stable structures (characteristic of the one-dimensional case) in the two-dimensional case is encountered only in as narrow class of the initial distributions. In most cases of random realizations of the initial inhomogeneity, the evolution of inhomogeneity proceeds rather smoothly.