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

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

Список литературы

1. Teunou, E. Batch and continuous fluid bed coating - review and state of the art./ E.Teunou, D. Poncelet.- Journal of Food Engineering, №53. - 2002. - P. 325-340.

2. Чуешов, В.И. и др. Промышленная технология лекарств: учебник в двух томах // Под ред. В.И. Чуешова- Х.: МТК-Книга, изд. НФАУ, 2002. -Т.2. -716 с.

3. Link, K.C. Fluidized bed spray granulation. Investigation of the coating process on a single sphere./ K.C. Link, E.-U Schlünder.- Chemical Engineering and Processing -№ 36. - 1997. - P. 443-457.

4. Меньшутина, Н. В. Моделирование и оптимизация тепло и массообмена на основе механики гетерогенных сред и неравновесной термодинамики в фонтанирующем слое: Дис... канд. техн. наук // РХТУ им. Д.И. Менделеева. - М.: РХТУ, 1985. - 174 с.

5. Корнеева А.Е. Моделирование атмосферной сублимационной сушки в аппаратах с активной гидродинамикой диссертация: Дис. канд. техн. наук // РХТУ им. Д.И. Менделеева. - М.: РХТУ, 2005. - 182 с.

УДК 519.711.3:541.12.012.3

П.А. Гуриков, А.В. Колнооченко

Российский химико-технологический университет им. Д.И. Менделеева, Москва, Россия

ТРЕХМЕРНЫЕ КЛЕТОЧНЫЕ АВТОМАТЫ ДЛЯ МОДЕЛИРОВАНИЯ СТРУКТУРООБРАЗОВАНИЯ ГЕЛЕЙ

Mathematical model based on three-dimensional stochastic cellular automata for simulation of the structural dynamics in gels were considered. The transition function and the criterion of the solid-liquid equilibrium in the model system were proposed. A software implementation of the model was developed. Identification of model parameters was implemented by genetic algorithm.

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

Развитие теоретических основ моделирования и новых подходов к моделированию сложных процессов, имеющих пороговых характер, является важным аспектом современных исследований. Клеточные автоматы позволяют моделировать сложное поведение динамической системы посредством разбиения ее на подсистемы, подчиняющиеся простым правилам [1].

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

Рассмотрим гетерогенную систему, в которой происходит образование твердой фазы, и выделим в ней кубический объем V. Рассечем этот объем N — 1 взаимно пер-

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

Всякая ячейка может находиться в одном из двух состояний: «раствор» или «кристалл». Концентрация компонента в растворе может меняться от 0 до плотности этого компонента в твердой фазе р. В состоянии «кристалл» концентрация совпадает с р. Далее для краткости будем говорить о «жидких» и «твердых» ячейках. Введем также понятие ближайших соседей, под которыми понимают ячейки, имеющие с данной общую грань. Если при определении координат ближайших соседей, координаты выходят за пределы [l, N] применяются т. н. периодические граничные условия:

i' = i mod( N +1)

j' = j mod( N +1) (1)

к' = к mod( N +1),

где ii, j', к' - новые координаты, i, j, к - старые координаты, mod - операция приведения по модулю (остаток от деления).

Массоперенос компонента между жидкими ячейками осуществляется за счет диффузии, описываемой законом Фика:

^ = D 0t

f Л

о c о c о c

+-2+^

ydx 0y 0z у

(2)

где Б - коэффициент диффузии, с - концентрация.

Начальные условия (состояние ячеек и концентрация компонента в жидких ячейках) и граничные (на границе раствор-кристалл массопереноса нет) считаются известными:

с 0(i, J, к) = J, к) (3)

Лс(г, J, к) = 0

Для численного решения уравнения диффузии используется явная разностная схема первого порядка аппроксимации совместно с условиями (1) и (3):

с"+1(/, J, к) = с" (/, J, к) + {с" (/ +1,;, к) + с" (/-1,J, к) +

Лх (4)

+ с" (/,^ +1,к)+с" (/, ^-1,к) + с" (/, J, к +1) + с" (/, J,к -1) - 6с" (/, к))

которая является устойчивой при выполнении условия

Л 1

-<--(5)

Лх2 6Б

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

Кристаллизация жидкой ячейки и растворение твердой в каждый момент времени Лt есть случайные события, происходящие с некоторыми вероятностями. Очевидно, что вероятность кристаллизации р жидкой ячейки с координатами (/, J, к) должна быть равна 0 при с(г, J, к) < С0 и достигать максимального значения, равного 1, при с(г, J, к) = р. Концентрация С0 есть концентрация компонента в растворе, в условиях

равновесия жидкой и твердой фаз. Вероятность растворения q твердой ячейки с координатами (г, 7, к), зависит от средней концентрации в ближайших соседних

с (г, 7, к) = 1 ^ с(г', 7', к') и равна 0 при с (г, 7, к) = р (все соседние ячейки твердые), а

при с (г, 7, к) < С0 равна 1. Наиболее простыми зависимостями р = р(с) и q = q(c) яв-

"0

ляются линейные:

0, с(г, 7, к) < С,

1, с (г, 7, к) < С(

0

р = {с(г, 7, к) - С0 q = {с (г, 7, к)-р (6)

Р-С0

c(i, 7;к) > С0

, С0 < с (г, j, к) <р С0 -Р

В момент реализации события «кристаллизация» концентрация вещества, как правило, меньше плотности кристалла, поэтому недостающее вещество должно поступить из соседних ячеек. Пусть кристаллизуется ячейка (г, 7, к) с концентрацией компонента с(г, 7, к). Тогда в эту ячейку должно поступить р — с(г, 7, к) вещества дополнительно. Рассмотрим всех жидких ближайших соседей ячейки (г, 7, к) - ячейки первого слоя. Вычислим суммарную концентрацию компонента ^ в них. Если < р — с(г, 7, к), то найдем суммарную концентрацию $2 в ячейках второго слоя (являющихся ближайшими соседями ячеек первого слоя) и т. д. до тех пор, пока для слоя с номером п выполнится условие ^ + £2 +... + £и > р — с(г, 7, к) . В этом случае в ячейках, принадлежащих слоям 1,2,..., п — 1, положим концентрацию компонента равной нулю, а в слое с номером п заберем вещество из ячеек пропорционально их изначальным концентрациям:

с(1\ 7, к') = с(г\ 7, к')8п - + +...+ ^п-1), (7)

^п

где с(г' , 7' , к' ), с(г', 7', к' ) - концентрация компонента до и после кристаллизации ячейки (г, 7, к) соответственно.

Если в ячейке реализуется событие «растворение», то изменяется лишь ее статус, а массоперенос в окружающий раствор происходит за счет диффузии по уравнению (4).

Кроме диффузии, растворения и кристаллизации на каждом шаге по времени At одиночные твердые ячейки способны перемещаться с вероятностью р вдоль одной из координатных осей влево или вправо - в позицию одного из жидких ближайших соседей (направление смещения и координатная ось выбирается случайным образом). Если смещение произошло за границу куба V, то применяются периодические граничные условия (1)

Состояние равновесия жидкой и твердой фаз характеризуется равенством скоростей образования и растворения твердой фазы, поэтому число твердых ячеек в состоянии равновесия должно испытывать случайные флуктуации около среднего значения (рис. 1). Этот факт используется для формулировки критерия окончания вычислительной процедуры. На протяжении некоторого числа итераций K (Х = 50 ^ 100) вычисляется среднее значение и дисперсия количества твердых ячеек (^ ,а2) . По прошествии следующих K итераций вычисляются новые оценки (а2 ,а2). Если считать, что число твердых ячеек распределено нормально, то задача сводится к проверке двухсторонней гипотезы о равенстве средних Н ,а1=а2. При условии, что равенство дисперсий и^ и

а рпоп не предполагается, гипотеза принимается на уровне значимости Р ( Р = 0,05 0,10 ), если

2 2 а: + а:

К

< г(п, Р),

(8)

где г(п, Р) - процентная точка распределения Стьюдента с числом степеней свободы, равном п :

п = -А

К

С _2 л

С -2 у

(9)

с а

+

К К

К —1

Рис. 1. Изменение числа твердых ячеек в процессе установления равновесия

Число п, вообще говоря, дробное. Процедура вычисления процентных точек с дробным числом степеней свободы приведена в [2].

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

Описанная выше математическая модель программно реализована. Численные значения параметров модели: размера куба N, плотности твердых ячеек р, коэффициента диффузии Б, равновесной концентрации С0, вероятности рт, шагов по времени Аг и координатам Ах = Ау = &, а также значения начальных концентраций в жидких ячейках ср(1,7, к) и координаты твердых ячеек читаются из входного файла. Правила перехода на каждой итерации даются приведенными выше формулами.

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

а — а2

2

щадь внутренней поверхности Sin, средний радиус пор г и т. д.) вычисляются после установления равновесия в модельной системе. Расчетные и экспериментальные зна-

exp oexp -exp

чения pg , S ^ , 'р предварительно нормируют.

Ясно, что вычисленные параметры аэрогелей будут зависеть от исходных параметров модели, в частности от коэффициента диффузии и начальной концентрации кремневой кислоты. Для определения таких значений D и p(i, j, к) = р = const, при которых экспериментальные значения минимально отличаются от вычисленных, применяется генетический алгоритм.

Все значения D и р предварительно нормируют на D0 и р. Значение D0 получают из условия (5). На первом шаге с помощью генератора псевдослучайных чисел формируют M уникальных пар (D1, р1), т. н. поколение родителей. Для каждой пары производят расчет, вычисляют параметры р , Sjn, Г и функцию пригодности R :

. R = ki(pgxp-Pg )2 + k2(S exp -Sin )2 + кз(Грехр - Гр )2 +... (10)

Коэффициенты значимости к, к2,... выбирают таким образом, чтобы ^ ki = 1. В таком

i

виде значения функции R изменяются от 0 (при точном соответствии параметров модели экспериментальным) до 1. Новое поколение (D2, р2) генерируют, отбирая случайно и независимо значения D и р из родительских пар, пропорционально их пригодности. С вероятностью pmMÍ (pmut ~ 0,05) произвольный параметр пары подвергают «мутации» -случайному изменению на псевдослучайное число из интервала [0, l]. После этого процедуру расчета и вычисления пригодности повторяют до тех пор, пока в одном из поколений

будет найдена пара (Dopt, popt) с пригодностью близкой к нулю.

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

Список литературы

1. Ванаг, В.К. Исследования пространственно распределенных динамических систем методами вероятностного клеточного автомата/ В.К. Ванаг.- Усп. физ. наук, 1999.- Т. 169.- № 5. -С. 481-505.

2. Большев, Л.Н. Таблицы математической статистики/ Л.Н. Большев, Н.В. Смирнов. -М.: Наука, 1983. - 416 с.

3. Hoshen, J. Percolation and Cluster Distribution. I. Cluster Multiple Labeling Technique and Critical Concentration Algorithm./ J.Hoshen, R.Kopelman.- Phys. Rev. B. 1976. V. 14. P. 3438-3445.

УДК 543: 615.2/.3

А.А. Войновский, А.И. Козлов

Российский химико-технологический университет им. Д.И. Менделеева, Москва, Россия

ИННОВАЦИОННАЯ ТЕХНОЛОГИЯ КОНТРОЛЯ КАЧЕСТВА ПРОИЗВОДСТВА ТВЕРДЫХ ЛЕКАРСТВЕННЫХ ФОРМ