CC BY
36
УДК 51-74
DOI:10.18413/2411-3808-2018-45-1-103-110
ПРИМЕНЕНИЕ МЕТОДА ДИСКРЕТНЫХ ЭЛЕМЕНТОВ ДЛЯ МОДЕЛИРОВАНИЯ
ВИБРОУДАРНЫХ СИСТЕМ
APPLICATION OF THE METHOD OF DISCRETE ELEMENTS FOR MODELING OF
VIBRO-IMPACT SYSTEMS
О.А. Верзилина O.A. Verzilina
Старооскольский технологический институт им. А.А. Угарова (филиал) федерального государственного автономного образовательного учреждения высшего образования «Национальный исследовательский технологический университет «МИСиС», Россия, 309516, г. Старый Оскол, микрорайон им. Макаренко, 42
Stary Oskol technological Institute. A. A. Ugarov (branch) of Federal state Autonomous-th educational institution of higher professional education "national research technological University "MISIS", 42 neighborhood to them. Makarenko, Stary Oskol, 309516, Russia
E-mail: [email protected]
Аннотация
Рассматривается задача применения метода дискретных элементов к задаче моделирования виброударных систем. Описаны основные этапы моделирования. Ввиду очень большой размерности системы обоснован и разработан метод моделирования с распараллеливанием вычислений при помощи технологии NVIDIA CUDA, который может быть полезен конструкторам авиационных деталей в связи со значительным сокращением времени моделирования. Проведено моделирование процесса виброударного упрочнения стыкового профиля нервюры крыла и получены значительные результаты сокращения времени, при этом сохраняется точность вычисления технологических параметров.
Abstract
In the current article considers the objective of applying the method of discrete elements to the problem of modeling of vibro-impact systems. The author describes the main stages of modeling. The modeling stage, which requires large computational operations and costs of machine time, is a modeling and numerical calculation of dynamic characteristics based on the discrete element method(DEM).Direct sequential calculations of these algorithms due to the very large dimensionality of the system (10,000 -40,000 equations) require significant machine time. In this regard, the method of modeling with parallelization of calculations using technology NVIDIA CUDA, which can be useful to designers of aircraft parts in connection with a significant reduction in simulation time, is justified and developed. The simulation of the process vibration shock hardening of the butt profile of the wing rib is carried out and significant results of time reduction are obtained, while the accuracy of the calculation of technological parameters is preserved.
Ключевые слова: метод дискретных элементов, виброударное упрочнение, NVIDIA CUDA, параллельны вычисления, сокращение времени.
Key words: discrete element method, vibro hardening, NVIDIA CUDA, parallel computing, reducing the time.
Введение
В предлагаемой статье представлена возможность компьютерной реализации моделирования работы установки поверхностного виброударного упрочнения. Использование этой технологии особенно важно для деталей, полученных литьём или
штамповкой по моделям с недоступным расположением внутренних поверхностей, эксплуатируемых в тяжёлых условиях: втулки винтов авиадвигателей, корпуса турбонасосных агрегатов и др. При этом, элемент конструкции помещается в контейнер, заполненный инструментальной средой из мелких абразивных частиц и промывочной жидкостью. Внешние вибраторы вызывают колебания контейнера, в результате циркуляционного перемещения рабочей среды и упруго-ударных контактов частиц с деталью, её поверхности придаются необходимые квазиупругие и вязкопластические свойства.
Основные этапы моделирования. Математическая модель
Компьютерное моделирование процесса виброупрочнения проведём методом дискретных элементов(МДЭ) [Белов, 2015; Williams, 1985; Cundall, 1979], сочетающим методы молекулярной динамики и теории упругопластических деформаций при ударе [Дорофеенко, 2008.]. Пространственное моделирование, предполагает подробное рассмотрение процесса в плоскостях, совпадающих с траекторией колебаний контейнера [Копылов, 2011] и последующей интерполяцией полученных результатов. В поперечных (продольных) сечениях рассчитываются координаты и скорости частиц инструментальной среды, определяются частота и энергии соударений частиц и обрабатываемой детали [Копылов, 2011]. Определяются съем металла, шероховатость, остаточные напряжения и т.д.
Выделим основные этапы моделирования работы установки. Этап 1. Построение трехмерной геометрической модели системы контейнер-деталь-инструментальная среда (рис.1).
Рис. 1 Трехмерная геометрическая модель Fig. 1. Three-dimensional geometric model
z
x
Этап 2. Определение набора характерных сечений в плоскостях, совпадающих с траекторией колебаний системы.
Этап 3. Аппроксимация границы. Контуры сечений (рис.2) и контейнера разбиваются на малые криволинейные отрезки, приближенно равные размеру частиц инструментальной среды. Каждый из них аппроксимируется линейной, квадратичной или более подходящей функцией - сплайном [Петряев, 2003; Петряев, Шевцов, 2001; Шевцов, 2001, и др.].
Рис. 2. Контуры поперечных сечений детали Fig. 2. The contours of the cross sections of the part
Этап 4. Моделирование загрузки контейнера. Определение исходного положения частиц инструментальной среды (рис.3).
Рис. 3. Заполнение контейнера с деталью(Б-Б) инструментальной средой Fig. 3. Filling the container with the part (B-b) of the tool environment
Этап 5. Моделирование и численный расчёт динамических характеристик на основе метода дискретных элементов(МДЭ) [Munjiza, 2004; Bicanic, 2004; Williams, 1999; Pande, 1990; Pande, 1990; Williams, 1989; Williams, 1985; Cundall, 1979].
Этап 6. Для каждого сплайна сечения детали рассчитываются динамические и технологические параметры: скорости и углы соударений, сила и продолжительность контакта, кинетическая энергия соударений, средняя высота микронеровностей (шероховатость), величина сжимающих остаточных напряжений [Копылов, 2011].
В завершении проводится визуализация трехмерного распределения динамических и технологических параметров процесса, анализ и статистическая оценка формирования шероховатости и остаточных напряжений, с выявлением проблемных участков и указанием их расположения.
Опишем подробнее центральный с теоретической и практической точки зрения, требующий больших вычислительных операций и затрат машинного времени, пятый этап. Дискретный материал инструментальной среды представляется совокупностью N упругих частиц (шариков) диаметра D i. Обозначим ri - радиус-вектор центра i -ой частицы, ñi -
её скорость, m¡ - массу, фг- - угол поворота вокруг центра, I¡ - момент инерции и шг- -
угловую скорость. Движение (поступательное и вращательное) каждой частицы полностью описывается основными уравнениями динамики [Белов, 2015; Копылов, 2011; Копылов, 1999]:
Г (t ) i = ñi (t )
U (t ) = — Fi (t ) mi
Ф i =®i(t )
i, j e [1, N],
(1)
где t - время,
® i = y Mi (t ) 1 i
Fii (t) = Z Fjn (A-, U-, t) + Z F4% (Aj-, ñj, t)
j: j * i, A ij > 0 j: j *i, Aj > 0
сумма сил,
действующих
M, (t ) = D S
на
контактах
1-й
частицы суммарный
с
соседними
7-ми,
Fij (Aij, ñij, t) -ч(t)
j: j *i, A j > 0
относительно центра i-й частицы. Здесь U- = U¡ — U- = (u
момент контактных
сил
ijn;Uij )
контактирующих частиц,
A = 1( D + D )—— j
) - встречная скорость величина их «перекрытия»,
<
(А;й;,, (Ау, й;, I) - нормальная и тангенциальная составляющие силы
контактного взаимодействия /-й и у-й частиц; т у - единичный вектор касательной в точке контакта.
Контактные силы Г ^) рассчитываются с учетом упруго-вязкого взаимодействия /й частицы со всеми соседними и представляют собой сумму сил отталкивания (2) и трения (3). В статье реализованы формы взаимодействия, использующее классический закон Кулона, контакт Уолта-Людинга, модели вязкого трения Фойхта-Кельвина, модели Финни-Рабиновича. При этом:
г ■ — г.. _ 3_
Щ =
г ■ — г-
3 г
г; = •(— к1 ■ Аз + к2 ■ щ]п)• пз, (2)
- единичный вектор нормальной компоненты, к1 и к2 - коэффициенты
модели нормального взаимодействия Герца-Кувабара-Коно
-т
Г; = —Ц
~п
й.
ут
ут
■Т; , (3)
ц - коэффициент, характеризующий свойства рабочей среды.
Границы контейнера и обрабатываемой детали разбиваются на малые участки, которые аппроксимируются известными функциями - сплайнами, и представляются совокупностью крупных псевдочастиц большой массы, что позволяет расширить систему (1), включив в контактные взаимодействия псевдочастицы. При этом, если М - число добавочных псевдочастиц, то в системе (1) ; е[1, N+М] и силы контактного
взаимодействия (А;, й;, ¿), Г;х (А;, й;, ¿) при ; > N +1 рассчитываются по формулам
несколько отличающимся от (2) и (3), учитывающих заданные внешние режимы движения контейнера и детали. Такой подход даёт возможность определять динамические и технологические параметры установки решая задачу Коши для системы (1). Методы численного решения подобных систем хорошо разработаны. Для программной реализации было выбрано комбинированное использование конечноразностного многошагового метода Адамса четвертого порядка и явного двухшагового. Общая схема подобного подхода для системы
Щ = Г (Щ, I) , Щ (Г0) = Щ
имеет следующий вид: предиктор
1/к + ГЧ 24 '(55Г (Щк, 1к ) — 59Г (Щк—1, 1к—1) + 37 г (Wk—2,1к—2) — 9Г (Щк—3,1к—3)) = Гк ;
корректор
= Жк
Г(Щк—1,%—1), Г(Щк_2, tk-2), Г(Щк—з,1к—3) - значения на предыдущих временных шагах, получаемые методом двухшагового интегрирования(4-5)
Щк+1/2 =Жк- Г(Щк ,'к ) ■ 42, (4)
Щ+1 = Жк + Г (Щк+12,'к+12) ■ 4. (5)
ЖкБ+1 = Жк + А/ 24 ■ (55Г Щ, гк ) — 59 Г Щ—1, гк—1) + 31Г Щ—Ъ гк—2) — 9Г Щ—3, гк—3)) = • Щы = Щ + V 24 ■ (Жк+1Б, гк+1) + 19Г (Жк, гк)—5Г (Щк ч, гкч) + Г (Щк—2, гк—2))= ТкК
Распараллеливание вычислений
Прямое последовательное проведение расчетов по этим алгоритмам ввиду очень большой размерности системы (10000 - 40000 уравнений) требует значительных затрат машинного времени. Кроме того, каждое передвижение частиц инструментальной среды может приводить к серьёзным изменениям правых частей системы (1), что заметно
усложняет вычислительный процесс. Поэтому на первый план выходит задача разумной и эффективной организации расчетов. Тем более, что оправдано с технологической точки зрения, естественно добиваться работы модели в режиме близком к реальному времени. Однако и сама структура матрицы системы (1) подсказывает необходимость и возможность распараллеливания вычислений. Одним из основных моментов при этом является разбиение рабочей области на достаточно малые ячейки, рассмотрение системы (1) в каждой из них и постоянный пересчет входящих элементов с изменением индексирования. Техническая реализация распараллеливания достигается за счет перераспределения вычислений по нескольким графическим мультипроцессорам.
Для распараллеливания вычислений автором настоящей работы выбрана технология CUDA [Боресков, 2012; Боресков, Харламов, 2010 и др]. CUDA - это архитектура параллельных вычислений от NVIDIA, позволяющая существенно увеличить вычислительную производительность благодаря использованию GPU (графических процессоров). GPU состоит из мультипроцессоров (streaming multiprocessor, MP), каждый из которых содержит несколько CUDA-ядер (CUDA core) и общую для них память.
Одну и ту же функцию - ядро (kernel) параллельно выполняют большое количество потоков-нитей (thread) [Боресков, 2012; Боресков, Харламов, 2010 и др]. Потоки группируются в блоки (thread blocks), при этом каждый блок исполняется на одном мультипроцессоре, а его потоки - на CUDA-ядрах данного мультипроцессора. Блоки объединяются в решетку/сетку блоков (grid). Ядро выполняется на решетке из блоков. Размер блока и размер решетки блоков задается при вызове ядра (рис. 4).
Grid
Ыаск(0,0) Ыоск(0,1) Ыоск(0,п-1)
Ыоск(1,0) Ыоск(1,1) Ыоск(1,п-1)
Ыоск(т-1,0) \Ь1оск(т-1,1) Ыоск(ш-1, и-1)
1 \ V
Block
tltread(0,0) 1hread(0,t-l)
threaa(k-l,0) ЛтаЛ(к-1,1-1)
Рис. 4. Иерархия нитей в CUDA Fig. 4. The hierarchy of threads in CUDA
В программе моделирования с распараллеливанием каждый этап интегрирования методом Адамса в момент времени t для численного решения задачи Коши реализован в виде отдельной функции-ядра, которая состоит из следующих ядер [Верзилина, 2016]:
-ядро copyBoundary - выполняет копирование параметров сплайнов (координат начала и конца, длины, нормали и т.д.), параметров псевдочастицы(позиции, скорости, угловые скорости, фазы), технологических параметров сплайнов(начала контакта, продолжительность контакта, скорость соударения, съем, наклеп, остаточные напряжения и т.д.) в момент времени time.
-ядро buildboundary преобразует координаты и скорости сплайнов и псевдочастиц согласно выбранному пользователем на этапе построения модели закону движения сплайновой границы.
-ядро predictor - ядро, реализующее схему предиктора метода Адамса для положений, скоростей, угловых скоростей и фаз ансамбля частиц инструментальной среды.
-ядро makePosVel - вспомогательная функция ядро, осуществляющее копирование позиций и скоростей частиц из двумерного массива в одномерный, для использования в дальнейшем в процедуре индексирования.
-ядро calcHash вычисляет в какую ячейку области разбиения попадает частица инструментальной среды на основе ее координаты.
-ядро sortParticles осуществляет сортировку частиц на основе информации о их номере и номере ячейки, в которую данная конкретная частица попала.
-ядро findCellStart поиск начала любой заданной ячейки в отсортированном списке. В результате выполнения ядер calcHash, sortParticles и findCellStart в дальнейшем при поиске взаимодействий только те частицы, которые непосредственно находятся в выбранной ячейке либо рядом с ней.
-функция D, выполняемая на GPU, вычисляет значение Fк, разделятся
дополнительно на 3 ядра:
-ядро D_Device - расчет взаимодействий между частицами инструментальной среды. -ядро D Spl Device - расчет взаимодействий между частицами инструментальной среды и сплайновой границей
-ядро Overrading - осуществляет пересчет позиций и скоростей и фаз с учетом полученных сил в момент времени +1.
-ядро corrector реализует схему корректора для положений, скоростей, угловых скоростей и фаз ансамбля частиц инструментальной среды.
-ядро swapArrays, выполняющее процедуру переопределения функций для следующего шага интегрирования.
Алгоритм вызова описанных выше ядер [Верзилина, 2016] представлен на рис. 5.
calcHashl)
|нет |
Рис. 5. Алгоритм моделирования шага интегрирования с распараллеливанием Fig. 5. The simulation algorithm of the step of integration with parallelization
Заключение
В результате распараллеливания вычислений процедуры интегрирования произведены замеры времени моделирования последовательной программы и программы моделирования, реализованной с распараллеливанием.
Исследование сокращения времени моделирования для тестовых моделей с числом частиц инструментальной среды от 500 до 10000 проводилось на видеокартах: Nvidia GeForce GT 640, Nvidia GeForce 740m, Nvidia GeForce 970. Максимально ускорение моделирования наблюдалось на видеокарте Nvidia GeForce 970 и на модели с числом частиц инструментальной среды 9000 составило 17 раз.
Используя программный комплекс с распараллеливанием вычислений, проведено моделирование процесса виброударного упрочнения стыкового профиля нервюры крыла ИЛ-96-300М содержащего 1795 сплайнов и 10000 частиц инструментальной среды, общий вид которой представлен на рис. 6.
Рис. 6. Общий вид стыкового профиля нервюры крыла Fig. 6. General view of the butt profile of the rib wing
Максимальное сокращение времени моделирования составляет 25 раз и позволяет получить параметры поверхностного слоя детали за 26 мин в параллельной реализации программы в сравнении с последовательным выполнением в течение 660 мин.
Список литературы References
1. Белов В.В., Образцов И.В., Иванов В.К., Коноплев Е.Н., 2015. Компьютерная реализация решения научно-технических и образовательных задач: учебное пособие. Тверь, ТвГТУ, 108.
Belov V.V., Obraztsov I.V., Ivanov V.K., Konoplev E.N., 2015. Computer implementation of solutions of scientific-technical and educational challenges: a tutorial. Tver, Tvgtu, 108.
2. Боресков А.В., Харламов А.А., 2010. Основы работы с технологией CUDA. М., ДМК Пресс, 232.
Boreskov A.V., Kharlamov A.A., 2010. Basics of CUDA technology. M., DMK Press, 232.
3. Боресков А.В. и др., 2012. Предисл.: Садовничий В.А. Параллельные вычисления на GPU. Архитектура и программная модель CUDA: Учеб. Пособие. Москва, 336.
Boreskov A.V. et al., 2012. Foreword.: Sadovnichy V.A.. Parallel computing on the GPU. Architecture and programming model CUDA: Proc. Benefit. Moscow, 336.
4. Верзилина О.А., 2016. Взаимодействие программных платформ DELPHI и C++ в рамках исследования двумерной виброударной многомассовой системы с распределенными параметрами с целью сокращения времени моделирования. Современные наукоемкие технологии. 10(2): 233238.
Verzilina O.A., 2016. The Interaction of software platforms DELPHI and C++ in the study, two-dimensional vibro-impact multi-mass system with distributed parameters with the aim of reducing time simulation. Modern high technologies. 10(2): 233-238.
5. Верзилина О.А., 2016. Разработка эффективного алгоритма исследования двумерной многомассовой системы с распределенными параметрами. Фундаментальные исследования. 11(2): 255-259.
Verzilina, O.A., 2016. Development of an efficient algorithm for the study of two-dimensional multi-mass system with distributed parameters. Fundamental research. 11(2): 255-259.
6. Дорофеенко С.О., 2008. Моделирование сыпучих сред методом дискретных элементов: дис. на соискание уч. степ. к-та физ.-мат. наук. Черноголовка: РАН Институт проблем химической физики, 114.
Dorofeenko S.O., 2008. Simulation of granular media using discrete element method: dis. on competition Uch. step. rand. Fiz.-Mat. sciences'. Chernogolovka: RAS Institute of problems of chemical physics, 114.
7. Копылов Ю.Р., 1999. Виброударное упрочнение. Воронеж, ВИМВД, 386.
Kopylov Y.R., 1999. Vibro hardening. Voronezh, WINGD, 386.
8. Копылов Ю.Р., 2011. Динамика процессов виброударного упрочнения: монография. Воронеж, ИПЦ «Научная книга», 568.
Kopylov Y.R., 2011. The Dynamics of vibroimpact hardening: monograph, Voronezh, CPI "Nauchnaya kniga", 568.
9. Петряев А.А., 2003. Компьютерное моделирование динамики гранулированных сред в вибрационных технологических машинах. Ростов-на-Дону, 218.
Petryaev A.A., 2003. Computer simulation of the dynamics of granular media in vibration technology machines. Rostov-on-don, 218.
10. Петряев А.А., Тамаркин М.А., Шевцов С.Н., 2001. Технологические приложения нестационарных задач динамики бинарных гранулированных сред. Труды междунар. конф. Изд. ДГТУ, Ростов-на-Дону, 1: 226-228.
Petryaev A.A., Tamarkin M.A., Shevtsov S.N., 2001. Technological applications of nonstationary problems of the dynamics of binary granular media. Proceedings of the international. Conf. Ed. DSTU, Rostov-on-don, 1: 226-228.
11. Петряев А.А., Тамаркин М.А., Шевцов С.Н., 2001 Прямое моделирование динамики гранулированных сред в некоторых технологических процессах. 8-й всероссийский съезд по теоретической и прикладной механике. Аннотации до-кладов. Изд. УрО РАН, Екатеринбург, 485.
Petryaev A.A., Tamarkin M.A., Shevtsov S.N., 2001 Direct modeling of the dynamics of granular media in certain technological processes. 8-th all-Russian Congress on theoretical and applied mechanics. Annotations up-to-treasure. Ed. UrO ran, Yekaterinburg, 485.
12. Шевцов С.Н., 2001. Моделирование динамики гранулированных сред при вибрационной отделочно-упрочняющей обработке. Ростов-на-Дону, 287.
Shevtsov S.N., 2001. Modeling the dynamics of granular media under vibration finishing and strengthening treatment. Rostov-on-don, 287.
13. Cundall P.A., O.D.L. Strack,1979. A distinct element model for granular assemblies. Geotechnique, 109.
14. Williams J.R., Hocking G., Mustoe G.G.W., 1985. The theoretical basis of the discrete element method. NUMETA.
15. Pande G., Beer G., Williams J.R., 1990. Numerical modeling in rock mechanics. John Wiley and Sons.
16. Williams J.R., Pentland A.P., 1989. Superquadric and modal dynamics for discrete elements in concurrent design. National science foundation sponsored 1st U.S.: Conference of discrete element methods. Golden, CO.
17. Williams J.R., Mustoe G.G.W., 1992. 2nd international conference on discrete element methods. IESL Press.
18. Williams J.R., O'Connor R., 1999. Discrete element simulation and the contact problem. Archives of computational methods in engineering. 279-304.
19. Bicanic N., 2004. Discrete element methods in Stein, de borst. Hughes encyclopedia of computational mechanics. 1.
20. Munjiza A., 2004. The combined finite-discrete element method.Wiley.
