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

УДК 51-74:519.711:519.714:666.972.7

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

А.Н. Бормотов1, И.А. Прошин1, А.В. Васильков2

Кафедра «Автоматизация и управление», ГОУ ВПО «Пензенская государственная технологическая академия» (1); aleks21618@yandex.ru;

ОАО НПП «Рубин» (2), г. Пенза

Представлена членом редколлегии профессором В.И. Коноваловым

Ключевые слова и фразы: композиционные материалы; математическое моделирование; многокритериальный синтез; оптимизация структуры и свойств; структурообразование; управление качеством.

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

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

Привлечение средств анализа требует записи решения системы дифференциальных уравнений в замкнутой форме, что позволяет исследовать эволюцию только двухчастичных моделей. Изучение процесса флокулообразования и установление статистических закономерностей, характерных для реальных дисперсных систем, возможно с применением методов численного анализа, так как число частиц в моделях таких систем должно быть достаточно велико (порядка 102...103).

Алгоритм численного анализа

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

miri -^ ( - vi) = -'VUi, / = 1,N , (1)

где т/ - масса /-й частицы; г, = (х,, у,, гг) - ее координаты; N - количество частиц наполнителя; к - коэффициент, определяемый диссипативными свойствами дисперсионной среды; V/ - скорость дисперсионной среды в точке г/; П/ - потенциал в точке г, (в общем случае зависящий от характеристик дисперсионной среды, а также от характеристик и взаимного расположения всех остальных частиц системы).

Представим потенциал в правой части системы (1) в виде суммы

и, = U.b + U,

l,g

N

+ ZUj, p j=1 j*i

где и/ ь - потенциал взаимодействия с границами; и, ё - гравитационный потенциал; Пу р - потенциал парного взаимодействия; N - число частиц.

Потенциал парного взаимодействия обычно выбирается в форме потенциала

Ми

или потенциала Морзе

U

m

= Z akr,j к=1

m

■У=!

Ur

' Р()= Xак ехР(ькТ/, + ск)! к=1

а также в форме их суммы или произведения с функцией Гаусса [1, 2].

Численные значения коэффициентов ак, Ьк, ск, входящих в выражения для потенциалов, устанавливаются на основании предварительной информации о количестве и координатах минимумов, соответствующих положениям равновесия. Частным случаем потенциала Ми является потенциал Леннарда-Джонса

u к )= u 0

V л12

r0

V r,j

- 2

i \6 Л [0

V r,j J

(2)

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

Для описания эволюции лиофильных дисперсных систем, между частицами которых действуют только силы отталкивания, в выражении (2) достаточно сохранить второе слагаемое. Если интерес представляет только взаимное расположение частиц (но не время, за которое достигается это расположение), то можно принять

к

UP(,j)= I-1 >

где к - постоянная.

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

и (т/,ь )= и о ТЦ,

Т/,Ь

r

где г ь - расстояние от поверхности частицы до граничной поверхности; По, го -

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

Полагая скорость частицы Vг- = Гг- за новую переменную, перепишем уравнение (1) в виде системы 6Ы обыкновенных дифференциальных уравнений

Г,- = V,- ;

i = 1, N .

(3)

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

взаимодействия

Fj = F ('и ) = Î u ('и )

У -, у у у^ У1] модуль силы взаимодействия с границей

д

области ¥,ь = — П(г ь), силу тяжести Рг- г = , и, в предположении о лами-

дг ’ ’*

нарном движении дисперсионной среды, силу вязкого трения Рг- е = 6пцЯ, (V - Г/ ),

где Я/ - радиус /-й частицы (рис. 1).

Тогда система (3) примет вид

V, = g + -

m,

N

X Fj + F,b + \e j=1

V и ^

которую удобней представить в форме

х(/) = f (х, t), (4)

где х = (г1, Г2,..., , VI, V2,..., VN) - радиус-вектор системы частиц в 6Ж-мерном

фазовом пространстве; t - время.

Пусть начальные условия заданы в виде

х(0) = хо. (5)

1

Рис. 1. Силы, действующие на частицу и подлежащие учету в процессе моделирования

Простейшим методом численного решения задачи Коши (4), (5) является метод Рунге-Кутта первого порядка (метод Эйлера), состоящий в замене функции £ изменяющейся на отрезке [т;т + к], значением этой функции в точке т

х(т + h) = х(т) + hf (x, т) = х(т) + S(x, т, h),

(6)

где S(x, т, h) = hf (x, т) - функция шага явного метода Эйлера.

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

На практике всегда используется апостериорная оценка точности. При этом простейшим адаптивным алгоритмом является метод сгущающихся сеток. В соответствии с данным алгоритмом переход от точки х(т) к точке х(т + h) временной сетки первоначально выполняется за один шаг h и затем за два шага h/2:

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

Использование формул (7), (8) в совокупности с методом Рунге-Кутта четвертого порядка приводит к необходимости выполнения одиннадцати вычислений правой части системы (4) на каждом шаге интегрирования.

Альтернативным алгоритмом контроля шага является алгоритм, основанный на методе вложенных форм. При использовании данного метода на каждом шаге интегрирования выполняется шесть вычислений правой части системы (4)

Полученные значения используются с двумя наборами весовых коэффициен-

Xh (т + h) = х(т)+S(x, т, h) ;

(7)

Д(т + h) = x h (т + h)-Xh (т + h). 2

(8)

(

\

(9)

i=1

Д(т + h) = x(t + h)-x* (т + h).

(12)

При использовании уравнений (9) - (12) предварительная корректировка шага интегрирования выполняется в соответствии с соотношением

й(т + h) = h(x) 5/ 5 , (13)

\|Д(х + h)

где 5 - коэффициент, характеризующий точность решения.

Если новое значение шага h оказалось больше предыдущего h (точность на текущем шаге завышена), то полученное значение фазы х(т + h) принимается, и в

качестве нового значения шага выбирается h*. Если новое значение шага оказа-

лось меньше (точность занижена по сравнению с требуемой), то значение х(т + h) отбрасывается, и расчет повторяется для меньшего шага h . В качестве нормы |Д| в уравнении (13) используется наибольший из модулей

|д| = тах{Дj}, i = 1,6N . (14)

І

Расчетная схема, представленная уравнениями (9) - (14), положена в основу программного обеспечения численного анализа.

Реализация алгоритма

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

Решение системы (4) при начальных условиях (5) возможно выполнить с использованием универсальных пакетов (MathCAD, MatLab, Maple и др.), однако для подобного подхода характерна сравнительно невысокая вычислительная эффективность. Кроме того, процесс реализации расчетной схемы из уравнений (9) - (14) в рамках универсальных пакетов по затратам времени сопоставим с реализацией ее в автономном ПО.

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

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

В частности, предположение о постоянстве ускорения а позволяет заменить решение системы (4) вычислением положения частиц в соответствии с кинематическим уравнением

ат2

х(т) = х(0)+у(0)хн--------------------------------. (15)

2

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

Таким образом, характерные ограничения существующих пакетов ПО численного анализа не позволяют использовать их для проведения численных экспериментов с использованием модели (4).

Для исследования эволюции дисперсных систем разработано не имеющее аналогов ПО. Программный продукт является автономным (не требует для работы других пакетов численного анализа) и реализован на стандартном языке ANSI C (ANSI X3.158-1989) для операционных систем Windows NT/2000, а также POSIX-совместимых вычислительных платформ.

Укрупненная блок-схема ПО приведена на рис. 2.

Исполнение ПО начинается с синтаксического анализа командной строки, в которой передается имя управляющей программы.

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

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

Управляющая программа

PSDL - интерпретатор

инкапсуляция физических характеристик системы частиц

инкапсуляция начальных условии

инкапсуляция характеристик границ и дисперсионной среды

инкапсуляция объектов мониторинга

.0

I

о

S

ш

d

ф

о

X

(Я

о

Û.

(N

СО

С

а

■&

d

ф

Расчетное задание

_________А_________

Планировщик

расчет бинарных взаимодействий

решение системы обыкновенных дифференциальных уравнений

расчет унарных взаимодействий

мониторинг результатов решения

Визуальное представление и статистические характеристики системы

о

I о

Рис. 2. Архитектура разработанного ПО численного анализа

В процессе решения при каждом вызове функции вычисления правой части системы (4) выполняется ряд операций.

N (Ж -1)

1. Нахождение —---------- значений сил парного взаимодействия

2

¥г] = -¥р = ^( -Кг -Ку ), г = 1, N-1, } = ,

гг}

где Еу - модуль силы парного взаимодействия; Гу - вектор, соединяющий центры г-й и у-й частиц; Кг, К у - радиусы г-й и у-й частиц соответственно.

2. Нахождение N значений силы тяжести

^ = т § >

где тг - масса г-й частицы.

3. Нахождение N значений силы вязкого трения

¥г,е = бпцК (у - & ),

где п - вязкость среды; Кг, & - радиус и скорость г-й частицы соответственно;

V - скорость дисперсионной среды.

4. Нахождение сил взаимодействия с границами.

Все граничные плоскости заданы уравнениями в нормальной форме

п к г - рк = 0, к = 1, К ,

где Пк - единичный вектор нормали к-й плоскости; г = (х, у, г) - радиус-вектор

текущей точки; рк - расстояние от начала координат до к-й плоскости; К - число граничных плоскостей.

Вектор силы, действующей со стороны граничной плоскости, вычисляется на основании соотношения

®г'к,Ь = пкЕгк,Ь (йгк - Кг ) = пкЕгк,Ь (пкГг - р - Кг ),

где Фк - расстояние от г'-й частицы до к-й плоскости; Кг - радиус г' -й частицы (рис. 3).

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

Для определения граничной сферы достаточно указать ее центр г и радиус

Яг (где I = 1,Ь , Ь - число граничных сфер).

Вектор силы, действующей со стороны граничной сферы, вычисляется на основании соотношения

F i,b - Т--------------------т Fil,b (dil - Ri )- у--------- Fa ь (R- - |r- - гг-1- Rj ),

lrl - ril lrl - ril

где йц - расстояние от I-й частицы до 1-й сферы (рис. 4).

Сумма сил тяжести, вязкого трения, парного взаимодействия и взаимодействия с граничными поверхностями, определяет силу, действующая на -ю частицу и изменяющую ее положение:

O

Рис. 3. Схема взаимодействия частицы с граничной плоскостью

Рис. 4. Схема взаимодействия частицы с граничной сферой

N

K

L

Fi - X Fj + mig + Fi,e + X Fik,b + X Fil,b '

j-1 k-1 l-1

j

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

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

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

Для этого целесообразно использовать следующие показатели:

где Я4! - расстояния от поверхности г-й частицы (г = 1, N, N - число частиц) до поверхности четырех ближайших частиц (для г-й частицы усреднение проводится по числу ближайших); Я4ау, ^4^ - оценки среднего и среднего квадратичного

отклонений расстояния от г-й частицы (г = 1, N ) до четырех ближайших к ней

соответственно; Nni - число частиц, расстояние до поверхности которых (от поверхности г-й частицы) не превышает заданного значения; Nnav, Nnstd - оценки среднего и среднего квадратичного отклонений числа частиц, расстояние до которых от данной не превышает заданного значения соответственно; Nsк - число

частиц, находящихся в k-й подобласти ( к = 1, K , K - число подобластей) выпуклой оболочки всех N частиц; Nsav, Nsstd - оценки среднего и среднего квадратичного отклонений числа частиц, находящихся в к-й подобласти (к = 1, K ) выпуклой оболочки всех частиц соответственно; K - число подобластей.

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

Численный эксперимент

Моделирование эволюции производилось для систем, объемная доля дисперсных частиц в которых варьировалась от 0,01 до 0,15 (для лиофильных систем до 0,37). В высоконаполненных дисперсных системах расстояния между частицами достаточно малы (существенно меньше размеров частиц). Поэтому эволюция таких систем протекает в стесненных условиях, что увеличивает характерное время процесса структурообразования до величин, сопоставимых с продолжительностью процесса твердения вяжущего. В качестве потенциалов парного взаимодействия были использованы:

1) потенциал Леннарда-Джонса (модель лиофобной системы);

2) потенциал Ми, представленный единственным слагаемым, характеризующим силы отталкивания (модель лиофильной системы);

3) сумма потенциала Леннарда-Джонса с функцией Гаусса (модель лиофиль-ной системы при наличии сольватных слоев).

Полученной статистической информации достаточно для построения адекватных моделей и установления закономерностей влияния рецептурно-технологических факторов на процесс структурообразования дисперсных систем [4].

Работа выполнена при поддержке гранта АВЦП «Развитие научного потенциала высшей школы (2009-2011 годы)» на тему «Математическое моделирование и многокритериальный синтез строительных материалов специального назначения», рег. № 2.1.2/11488.

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

1. Мелькер, А.И. Самоорганизация и образование геликоидальных структур полимеров / А.И. Мелькер, Т.В. Воробьева // Физика твердого тела. - 1997. -Т. 39, № 10. - С. 1883-1889.

2. Моделирование процессов структурообразования дисперсных систем / А.Н. Бормотов [и др.] // Тр. Междунар. конф. «Идентификация систем и задачи управления» SICPRO’05 / Ин-т проблем упр. - М., 2004. - С. 700-724.

3. Numerical Recipes in C: The Art of Scientific Computing / William H. Press [at al.]. - Cambridge : Cambridge University Press, 1992.

4. Бормотов, А.Н. Имитационное моделирование деструкции и метод прогнозирования стойкости композиционных материалов / А. Н. Бормотов, И. А. Прошин, Е.В. Королев // Вестн. Ижев. гос. техн. ун-та. - 2010. - № 4. - С. 113-118.

Theoretical Bases of Computer Modeling of Disperse Systems Structurization

A.N. Bormotov1, LA. Proshin1, A.V. Vasilkov2

Department «Automation and Control»,

Penza State Technological Academy (1); aleks21618@yandex.ru;

Open Joint-Stock Company «Research-and-Production Association «Rubin», Penza (2)

Key words and phrases: composite materials; mathematical modeling; multicriteria synthesis; optimization of structure and properties; quality management; structurization.

Abstract: The paper presents the analytical solutions of the systems of interacting elements of the structure of composite materials, which are the basis of complex software simulation of structure formation. By means of computer simulation the evolution of disperse systems is studied, adequate mathematical models are built; and laws of influence of recipe-technology factors on the processes of structurization of disperse systems are established.

Theoretische Grundlagen der Computermodellierung der Strukturbildung der Dispersensysteme

Zusammenfassung: Es werden die analytischen Lösungen der Systeme der zusammenwirkenden strukturbildenden Elemente der Kompositstoffen, die die Grundlage des Komplekxes der Programmen der Computermodellierung der Strukturbildung sind, angeführt. Mit Hilfe der Mittel der Computermodellierung wird die Evolution der Dispersensysteme untersucht, es werden die adequaten mathematischen Modelle gebaut und es werden die Gesetzmäßigkeiten der Einwirkung der rezeptur-technologischen Faktoren auf die Prozesse der Strukturbildung der Dispersensysteme festgestellt.

Bases théoriques du modélage informatique des systèmes dispercifs formant des structures

Résumé: Sont citées les solutions analytiques des systèmes des éléments des matériaux composites coopérants et formant des structures. Par les moyens du modélage informatique est étudiée l’évolution des systhèmes dispercifs, sont construits des modèles mathématiques adéquats et sont établies les régularités de l’influence des facteurs réceptifs et technologiques sur les processus de la formation des structures des systèmes dispercifs.

Авторы: Бормотов Алексей Николаевич - кандидат технических наук, доцент, профессор кафедры «Автоматизация и управление»; Прошин Иван Александрович - доктор технических наук, профессор, заведующий кафедрой «Автоматизация и управление», ГОУ ВПО «Пензенская государственная технологическая академия»; Васильков Александр Васильевич - заместитель начальника отделения, ОАО «Научно-производственное предприятие «Рубин», г. Пенза.

Рецензент: Коновалов Владимир Викторович - доктор технических наук, профессор кафедры «Механизация животноводства», ФГОУ «Пензенская государственная сельскохозяйственная академия», г. Пенза.