CC BY
33
10000
8000
6000
4000
2000
1591122233
Рис. 5. Давление в поре без учета внутреннего
напряженного состояния:
- нет внешнего напряжения; ■
- растяжение
2000000 1600000 1200000 800000 400000 0
1 5 9 13 17 21 25 29 33 37
Рис. 6. Давление в поре: без учета внутреннего давления в механической задаче;----с учетом давления
P
0
t, с
На основании изложенного можно сделать следующие выводы. Близость графиков на рисунке 3 показывает хорошую согласованность результатов расчетов с аналитическим решением. Отклонение результатов расчетов от аналитического решения вблизи начального момента времени обусловлено тем, что принятый в качестве аналитического решения ряд дает на этом интервале плохое приближение [5]. Как видно из графиков на рисунках 4 и 5, учет напряженного состояния для двухмерной задачи диффузии необходим, причем под действием внешнего растягивающего напряжения диффузионный процесс существенно
ускоряется. График на рисунке 6 характеризует влияние диффузионного процесса на механическую задачу в околопоровой зоне. Из рисунков видно, что характер роста давления в поре для решения однородной задачи диффузии и для решения связной диффузионно-деформационной задачи различен. Таким образом, для анализа прочности плоских деталей металлических конструкций с дефектами необходимо учитывать как внутреннее напряжение в дефектах вследствие процесса диффузии, так и внешнее механическое напряжение вследствие процесса деформации. Представленный в работе алгоритм и программа учитывают оба этих фактора.
Литература
1. Андрейкив А.Е., Панасюк В.В., Поляков Л.И., Харин В.С. Механика водородного охрупчивания металлов и расчет элементов конструкций на прочность. Львов: препр. № 133. Физ.-мех. ин-т им. Г.В. Карпенко. 1987. 50 с.
2. Федотов В.П., Спевак Л.Ф. Модифицированный метод граничных элементов в задачах механики, теплопроводности и диффузии. Екатеринбург: УрО РАН, 2009. 164 с.
3. Нефедова О.А., Федотов В.П. Применение модифицированного метода граничных элементов для решения задач теплопроводности / Актуальные проблемы современной науки: тр. 5-го Междунар. форума молод. ученых. Самара: СамГТУ, 2010. С. 143-147.
4. Бреббия К., Теллес Ж., Вроубел Л. Метод граничных элементов. М.: Мир, 1987. 526 с.
5. Лыков А.В. Теория теплопроводности. М.: Высш. шк. 1967. 600 с.
References
1. Andreykiv A.E., Panasyuk V.V., Polyakov L.I., Ha-rin V.S., Lvov, Fiziko-mekhanicheskiy Institut imeni G.V. Karpen-ko, 1987, 50 p.
2. Fedotov V.P., Spevak L.F., Ekaterinburg, Uralskoe otdele-nie Rossiyskoy akademii nauk, 2009, 164 p.
3. Nefedova O.A., Fedotov V.P., 5th Mezhdunar. forum molod. uchenykh (5th Int. Forum for Young Scientists), Samara, 2010, pp. 143-147.
4. Brebbiya K., Telles Z., Wroubel L., Metod granichnykh elementov (Boundary element method), Moscow, 1987, 526 p.
5. Lykov A.V., Teoriya teploprovodnosti (The heat conduction theory), Moscow, 1967, 600 p.
УДК 66.071.6: 004.942-021
РАЗРАБОТКА АЛГОРИТМА ПАРАЛЛЕЛЬНЫХ ВЫЧИСЛЕНИЙ ДЛЯ ОПИСАНИЯ МАССОПЕРЕНОСА В ПОРЕ
(Работа выполнена в рамках ГК с Минобрнауки РФ №№ 16.740.11.0040, 11.519.11.4004 и гранта РФФИ № 11-08-91159 ГФЕН а)
А.Д. Поветкин; Чан Хыу Куе; Э.М. Кольцова, д.т.н.
(Российский химико-технологический университет им.. Д.И. Менделеева, г. Москва,
que_kola@m-ail.ru, lcGlts@muctr.Tu)
Рассматривается модель на основе метода молекулярной динамики, описывающая перенос газа в порах малого размера. Учитываются два типа переноса вещества, происходящего внутри поры: диффузия Кнудсена и молекулярная диффузия. Для описания движения и взаимодействия молекул рассматриваются два варианта. В первом молекулы движутся по законам классической механики и взаимодействуют друг с другом, а также со стенкой по принципу
абсолютно упругого удара. Принимается, что частицы могут сталкиваться со стенкой поры разными способами -зеркальным и диффузным, причем для каждого конкретного соударения способ определяется случайным образом, и соотношение количества соударений по обоим способам является одним из параметров модели. Во втором варианте взаимодействие молекул описывается с применением межмолекулярного потенциала взаимодействия Леннарда-Джонса. Расчет движения молекул проводится с привлечением алгоритмов параллельного вычисления. Организация параллельных вычислений осуществляется с применением технологии СЦЭА. Для хранения параметров частиц используется трехмерный массив. Число ячеек массива соответствует числу частиц в системе. В каждой ячейке массива хранятся параметры какой-либо частицы (координаты, вектор скорости, вектор ускорения и тип). Отдельно существует одномерный массив типов. Количество элементов массива равно количеству типов частиц, присутствующих в системе. В ячейках хранятся параметры каждого из веществ, присутствующих в системе, - масса и радиус частицы, коэффициенты потенциала межмолекулярного взаимодействия. Для упорядочения элементов в массиве применяется алгоритм сортировки пузырьком, приспособленный для параллельной реализации и расширенный для трехмерного случая. Приведены основные блок-схемы расчета.
Ключевые слова: молекулярная динамика, диффузия газа, диффузия Кнудсена, коэффициент диффузии, параллельное вычисление.
PARALLEL DESIGN ALGORITHM FOR DESCRIPTION OF THE MASS TRANSFER IN A PORE Povetkin A.D.; Tran Huu Que; Koltsova E.M., Ph.D.
(D. Mendeleev University of Chemical Technology of Russia, Moscow, que_kola@mail.ru, kolts@muctr.ru) Abstract. Here is described a model using molecular dynamics approach that describes gas transfer in small pores. Two types of the mass transfer that happen in a pore were taken into account: Knudsen diffusion and molecular diffusion. Two cases of interaction are mentioned that describe molecular movement and interaction. First case describes the variant, where molecules move according to classic mechanics law and where they interact with each other and a wall, on the principle of perfectly elastic collision. It is assumed that the particles can collide with the pore's wall in two ways: mirror and diffusion, and for each specific collision this way is defined at random fashion, and the ratio of collisions according to both ways presents one of the features of the model. The second case describes molecular interaction using intermolecular potential of Lennard-Jones interaction. Calculation of molecular movement was made with parallel algorithms. Parallel computing is made with the use of CUDA technology. Information about particle features is stored in 3-D array. Number of array cells corresponds to number of particles in the system. Each array cell contains information about particle: coordinates, velocity vector, speed up vector and its type. One dimension array of particle type exists separately. Number of array elements equals to number of types of existing particles. Cells contain information about each substance that is present in the system: particle's mass and radius, coefficient of molecular interaction potential. The array was arranged with bubble sort algorithm adapted for parallel operation and extended for 3-D case. The article contains basic computation flow charts.
Keywords: molecular dynamics, gas diffusion, Knudsen diffusion, diffusion coefficient, parallel computing.
Математическая модель на основе метода молекулярной динамики описывает перенос вещества в порах диаметром от 1 до 1000 нм, длиной от 0,01 до 200 мкм, при давлениях от 0,001 до 2 атм. При малых размерах пор перенос вещества [1] происходит как за счет столкновения молекул друг с другом (молекулярная диффузия), так и за счет соударений молекул со стенкой поры (диффузия Кнудсена), причем второй тип преобладает, так как вероятность соударения молекул друг с другом при низком давлении довольно мала. После задания начальных условий (длина, диаметр поры, температура, давление, молекулярная масса газа) производится расчет количества молекул в поре в соответствии с законом Менделеева-Кла-ЫАрУ
пейрона N =■
RT
где Na - число Авогадро; V -
объем поры.
Средняя скорость молекул газа рассчитывается из распределения Максвелла-Больцмана. После вычисления средней скорости молекул газа каждой молекуле присваивается определенная скорость с учетом вероятности из распределения Максвелла.
Взаимодействие молекул описывается с применением межмолекулярного потенциала взаимодействия Леннарда-Джонса
U {г) = 4в
где г - расстояние между центрами частиц; е -глубина потенциальной ямы; ст - расстояние, на котором энергия взаимодействия становится равной нулю. Параметры е и ст являются характеристиками молекул соответствующего вещества.
С использованием генератора случайных чисел и с учетом распределения давления по поре всем частицам (молекулам) присваиваются начальные координаты их местоположения. Далее координаты положения частиц меняются в соответствии с уравнениями движения.
При столкновении со стенкой пересчитывают-ся скорости с учетом зеркального и диффузного механизмов рассеяния. При зеркальном рассеянии частиц происходит абсолютно упругий удар, и частица отлетает от стенки под тем же углом, под которым она приближалась к ней. При диффузном рассеянии новое направление движения частицы определяется с помощью генератора случайных чисел. В обоих случаях модуль вектора скорости остается постоянным. Проверка частиц на столкновение проводится только для наиболее близко расположенных друг к другу частиц, для которых возникает наибольшая вероятность соударения.
6
Организация параллельных вычислений осуществляется с применением технологии CUDA. Преимущество использования технологии параллельного программирования в том, что при моделировании необходимо проводить множество однотипных операций над всеми частицами (проверка на столкновения, численное решение уравнений движения и др.). Для такой задачи хорошо подходят решения на основе принципа SIMD (single instruction-stream, multiple data-stream -один поток команд с несколькими потоками данных), когда вычислительное устройство выполняет операции сразу над несколькими элементами данных за один такт [2]. На этом принципе построена и используемая в данной реализации технология CUDA.
Для хранения параметров частиц используется трехмерный массив. Число ячеек массива соответствует числу частиц в системе. В каждой ячейке массива хранятся параметры какой-либо частицы: координаты, вектор скорости, вектор ускорения и тип. Отдельно существует одномерный массив типов. Количество элементов массива равно количеству типов частиц, присутствующих в системе. В ячейках хранятся параметры каждого из веществ, имеющихся в системе, - масса и радиус частицы, коэффициенты потенциала межмолекулярного взаимодействия. В ходе расчета при необходимости получения доступа к этим параметрам они берутся из ячейки с номером, соответствующим типу рассчитываемой частицы.
Алгоритм для сортировки частиц требует упорядоченности массива по координатам частиц. Для упорядочения элементов применяется простой алгоритм сортировки пузырьком [3], приспособленный для параллельной реализации, и расширенный - для трехмерного случая. Пусть имеется одномерный неупорядоченный массив a из n элементов. Алгоритм состоит в повторяющихся проходах по сортируемому массиву. За каждый проход элементы последовательно сравниваются попарно и, если порядок в паре неверный, выполняется их обмен (рис. 1). Проходы по массиву повторяются до тех пор, пока не окажется, что обме-
1 3 2 7 4 5 6 8
1 [Т
2 3 7
3 7 4
^7 5
^Г7 6 ^Р^ 8
7 8
1 2 3 4 5 6 7 8
Рис. 1. Работа алгоритма сортировки пузырьком
ны больше не нужны, это означает, что массив отсортирован.
При параллельной реализации каждый проход по массиву разделяется на 2 этапа. На первом попарно сравниваются элементы с номерами 2i и 2/+1 (нулевой с первым, второй с третьим и т.д.). На втором этапе происходит то же самое, но со сдвигом на 1 позицию, то есть сравниваются 2/+1 и 2/'+2 элементы (рис. 2).
Исходный массив
т- с га
О н О
X СМ
о
с^ с
1= га н
О
CN с
га н
О X О
о СМ
с^
1= с га
н
О
1 3 2 7 4 5 6 8
1 3 ||2 7 |4 5 ||6 8
1 3 2 7 4 5 6 8
1 8
3 2 | |7 4 || 5 6
1 2 3 4 1 5 6 8
1 2 ||3 4| |7 5 ||6 8
1 2 3 4 5 7 6 8
12 3 4 5 7 6 8
1 2 3 4 5 г-< 6 7 8
Рис. 2. Параллельная реализация пузырьковой сортировки
В результате после прохода массив становится частично отсортированным, как и в оригинальном алгоритме, однако в пределах одного этапа все сравнения могут выполняться параллельно. Понадобится п/2 потоков, каждый из них обработает одну пару элементов. Так как не будет происходить перекрывание данных, алгоритм не требует блокировок и все потоки могут выполняться одновременно.
Рассмотрим трехмерный массив с размерностью LxMxN. Причем, в отличие от одномерного случая, когда в ячейках массива хранилось всего одно число, в трехмерном случае каждая ячейка будет представлять собой вектор из трех координат у, z). Шаг сортировки для такого массива с применением вышеописанного одномерного алгоритма выбирается следующим образом. Очевидно, что трехмерный массив LxMxN можно представить в виде NxM одномерных массивов из L элементов. Для каждого из этих массивов следует выполнить проход алгоритма сортировки, причем для сравнения значения x брать из ячеек. После этого трехмерный массив будет частично отсортирован по одной из размерностей. Аналогично следует выполнить проход по двум другим размерностям, сравнивая соответствующие координаты, то есть для M-массивов берутся координаты у, а для ^массивов - координаты z. В итоге после такого прохода массив будет частично отсортирован по всем трем координатам.
В общем случае для полной сортировки массива проходы должны выполняться до тех пор, пока
не прекратятся перестановки элементов, и для изначально неупорядоченного массива таких проходов может быть довольно много. Однако для массивов, в которых значения уже частично упорядочены, особенно для таких, которые образованы несколькими перестановками элементов упорядоченной последовательности, алгоритм показывает хорошую производительность, так как количество проходов в таких случаях довольно мало.
Массив частиц следует сортировать после каждого расчета нового положения частиц в пространстве. При таком использовании алгоритма сортировки расположение частиц в ячейках массива будет приближенно соответствовать их положению в пространстве внутри поры в процессе моделирования, что является важным условием для корректного расчета взаимодействия частиц.
Проверка на столкновение выполняется только для тех частиц в системе, которые наиболее близко расположены друг к другу, и у них, соответственно, будет наибольшая вероятность соударения. Проверка частиц на соударение проводится в два этапа. На первом этапе массив разбивается на кубические блоки размером РхРхР (рис. 3). В каждом из этих блоков происходит полный перебор всех сочетаний из двух частиц, и для каждого выполняется проверка на соударение с последующим пересчетом скоростей при положительном результате проверки.
Второй этап аналогичен первому, но разбивка на блоки происходит со сдвигом, чтобы учесть столкновения, которые могли бы произойти между частицами, находящимися рядом, но попавшими на первом этапе в разные блоки. Сдвиг идет по всем трем размерностям массива на Р/2 позиций.
Для организации проверки частиц на соударение со стенкой для каждой частицы в массиве проверяется, выходят ли ее координаты за пределы поры. При положительном результате скорости частицы пересчитываются. При соударении, помимо пересчета скорости, сам факт регистрируется счетчиком соударений. Регистрируются длина свободного пробега частицы от последнего соударения и время, которое прошло с последнего соударения, после чего данные о времени и длине пробега обнуляются. Эти данные необходимы для расчета коэффициента диффузии D частиц в поре. Коэффициент диффузии рассчитывается по фор-
(v)(l)
муле D = , где (v) - средняя скорость частиц 6
в системе; (I) - средняя длина свободного пробега.
С целью проверки частиц на выход из поры и организации входа частиц в пору для каждой частицы в массиве проверяется, не превышает ли координата х частицы длину поры I. Если превышение есть, частица считается вылетевшей. Факт выхода из поры регистрируется счетчиком, а ячейка, которую заполняла вылетевшая частица, заполняется новой частицей. Ей присваиваются тип в соответствии с составом, координаты и скорости, как и при начальном заполнении, за исключением того, что координаты х, у новой частицы равны 0, так как частица находится на входе в пору.
Сила межмолекулярного взаимодействия рассчитывается только для близко расположенных друг к другу частиц. Так как расположение частиц в ячейках отсортированного массива примерно соответствует их положению в пространстве, для каждой ячейки массива в соседних ячейках, вероятнее всего, будут находиться частицы, наиболее близкие к частице в данной ячейке в пространстве поры. Для расчета сил межмолекулярного взаимодействия для каждой частицы в массиве берется совокупность соседей, находящихся в ячейках не более чем в R позициях от данной (рис. 4). Из рисунка 4 видно, что области частиц A и B перекры-
Рис. 4. Области вокруг частиц, для которых проводится расчет сил межмолекулярного взаимодействия
Рис. 3. Разбиение массива на блоки: А - без сдвига, Б - со сдвигом
ваются, при этом каждая из частиц попадает в область другой.
Рассчитывается только модуль силы, действующий на частицу в данной ячейке; для всех остальных ячеек выполняется такой же расчет в своей окрестности. Так как окрестности близко расположенных частиц перекрываются, для всех подобных пар будет корректно рассчитана сила
парного взаимодействия (каждая из частиц рассчитает свою силу).
Выходные статистические данные можно получить непосредственно (значения счетчиков соударений, счетчиков вышедших частиц), посчитав среднее арифметическое значение нужной величины для всех частиц в массиве (средняя скорость, средняя длина свободного пробега), а также вычислить по формуле (коэффициент диффузии).
Задание температуры, давления, радиуса и длины
поры, времени и шага расчета, объемных долей каждого вещества в смеси и параметров веществ
Сортировка частиц в массиве по положению в пространстве
I
Проверка частиц на столкновение друг с другом
I
Проверка частиц на столкновение со стенкой
Перемещение частиц в пространстве в соответствии с их скоростями
Ж
Обсчет входа и выхода частиц в пору через ее концы
Сбор статистики о частицах
(скорости, столкновения, вылетевшие частицы), расчет выходных параметров (выходной состав, коэффициент диффузии)
I
Вывод статистики на экран и в файл
т
Приращение текущего времени на 1 шаг
Рис. 5. Общая схема расчета
(Конец)
Рис. 6. Блок-схема сортировки
Все данные выводятся в удобной для пользователя форме на экран монитора или на запоминающее устройство.
Общий алгоритм расчета осуществляется следующим образом. После ввода пользователем параметров моделирования происходят определение размера массива (из уравнения Менделеева-Клай-перона) и начальное заполнение массива частицами. Частицам присваиваются тип в соответствии с заданным составом смеси (так как состав задается мольными долями, которые могут быть пересчитаны непосредственно в количество молекул, в итоге количество частиц разных типов в массиве отражает состав смеси), начальные пространственные координаты, начальные скорости (в соответствии с распределением Максвелла). Задаются шаг моделирования по времени и общее время моделирования, текущее время моделирования устанавливается на 0. После этого запускается процесс моделирования массопереноса в поре, который состоит в выполнении над массивом частиц повторяющейся последовательности операций (рис. 5):
- частичная сортировка массива одним проходом сортировки (рис. 6), в результате которой будут пересчитаны скорости для столкнувшихся частиц;
- фактическое выполнение моделирования (если проверка частиц на столкновение друг с другом не проводится) с учетом только кнудсе-новской диффузии, так как частицы не взаимодействуют друг с другом;
- проверка частиц на столкновение со стенкой, в результате которой будут пересчитаны скорости столкнувшихся со стенкой частиц;
- расчет сил межмолекулярного взаимодействия, в результате которого всем частицам будут приданы соответствующие ускорения; если расчет
не проводится, моделирование проходит без учета сил межмолекулярного взаимодействия;
- проверка частиц на выход из поры и организация входа частиц в пору;
- расчет нового положения частиц в соответствии с уравнениями движения, в результате чего частицам будет присвоено новое положение в пространстве;
- сбор статистики для определения выходных данных моделирования;
- увеличение текущего времени моделирования на величину шага и проверка, не превысило ли оно заданное время моделирования; при отрицательном результате повторение последовательности (то есть новый шаг по времени), при положительном расчет заканчивается.
Применение технологии организации параллельного вычисления CUDA для описания математической модели позволяет существенно быстрее производить расчеты для большого количества молекул N, которое в зависимости от условий (давление, температура, длина, диаметр поры) может варьироваться от 1 до 100 тыс. молекул по сравнению с однопоточной версией программы.
Литература
1. Albo S.E., Broadbelt L.J., Snurr R.Q. // AIChE Journal,
2006, Vol. 52, no. 11, pp. 3679-3687.
2. Таненбаум Э. Архитектура компьютера. СПб: Питер,
2007. 5-е изд. 844 с.
3. Левитин А.В. Алгоритмы: введение в разработку и анализ. М.: Изд-во «Вильямс», 2006. 576 с.
References
1. Albo S.E., Broadbelt L.J., Snurr R.Q., AIChE Journ., 2006, Vol. 52, no. 11, pp. 3679-3687.
2. Tanenbaum E., Arkhitektura kompyutera (Computer architecture), St. Petersburg, 2007, 844 p.
3. Levitin A.V., Algoritmy: vvedenie v razrabotku i analiz (Algorithms: An Introduction to Design and Analysis), Moscow, 2006, 576 p.
УДК 574.6.663.1
АЛГОРИТМЫ ОПТИМИЗАЦИИ НЕПРЕРЫВНОГО ПРОЦЕССА БИОСИНТЕЗА МОЛОЧНОЙ КИСЛОТЫ
Ю.Л. Гордеева, к.т.н.; Ю.А. Ивашкин,, д.т.н. (Московский государственный университет пищевых производств, l.s.gordeev@yandex.ru);
Л.С. Гордеев, д.т.н.
(Российский химико-технологический университет им.. Д.И. Менделеева, l.s.gordeev@yandex.ru)
Получены соотношения для расчета показателей непрерывного процесса биосинтеза молочной кислоты. В основе соотношений лежит математическая модель непрерывного процесса синтеза в ферментёре с перемешиванием. Особенность модели заключается в том, что для каждого компонента (субстрата, биомассы и продукта) записывается свое выражение для удельной скорости. В качестве критерия оптимальности используется продуктивность Qp по целевому продукту (молочной кислоте). При решении оптимальной задачи сначала оцениваются скорость протока Б и концентрация субстрата в выходном потоке, а затем рассчитывается концентрация субстрата на входе в аппарат. Эти показатели обеспечивают максимум продуктивности. Полученные соотношения использованы для разработки алгоритмов оптимизации непрерывного процесса биосинтеза молочной кислоты. Рассмотрены три варианта поста-
