Содержание к диссертации
Введение
1 Анализ исследований и практики моделирования многомассовых виброударных систем. Цель и задачи исследований 12
1.1 Развитие теории соударений 12
1.2 Анализ исследований многомассовых виброударных систем 15
1.3 Обзор методов моделирования многомассовых виброударных систем 19
1.4 Анализ существующих комплексов и программ моделирования многомассовых виброударных систем 24
1.5 Обзор программно-аппаратных средств для распараллеливания вычислений 27
1.6 Постановка целей и задач исследования 40
1.7 Выводы по первой главе 41
2 Описание объекта исследования. этапы моделирования двумерных виброударных систем. Структурная схема численного решения 42
2.1 Характеристика объекта исследования. Методы моделирования. 42
2.2 Описание основных этапов моделирования двумерных виброударных систем 49
2.3 Структурная схема численного решения модели 70
2.4 Выводы по второй главе 76
3 Определение резервов сокращения времени. разработка нового алгоритма вычисления с распараллеливанием при помощи технологии NVIDIA CUDA 77
3.1 Выбор тестовой модели многомассовой дискретной виброударной системы для определения затрат времени на вычислительные процедуры 77
3.1.1 Влияние количества дискретных частиц виброударной системы на затраты времени моделирования по ранее разработанным моделям и программам 78
3.1.2 Влияние законов ударного трения на затраты времени моделирования по ранее разработанным моделям и программам 82
3.2 Детализация функций, входящих в процедуру интегрирования и возможности их распараллеливания 83
3.2.1 Индексация области моделирования 83
3.2.2 Функция вычисления контактных сил - «D» 84
3.2.3 Функции Buildboundary, предиктор, корректор 86
3.3 Выбор видеокарты 87
3.3.1 Видеокарта Nvidia GeForce GT 640 87
3.3.2 Видеокарта Nvidia GeForce 740m 89
3.3.3 Видеокарта Nvidia GeForce 970 89
3.3.4 Сравнение GPU характеристик, выбранных видеокарт 92
3.4 Разработка нового алгоритма моделирования виброударных систем с распараллеливанием на основе программной платформы NVIDIA CUDA 95
3.4.1 Алгоритм компьютерного моделирования процедуры интегрирования с распараллеливанием 95
3.4.2 Алгоритм взаимодействия программных систем Delphi и C++ 102
3.5 Выводы по третьей главе 107
4 Разработка программы для реализации эффективного вычислительного метода моделировния с распараллеливанием на видеокартах NVIDIA 109
4.1 Выбор среды разработки 109
4.2 Взаимодействие Delphi и С++ 114
4.3. Описание комплекса программ с распараллеливанием вычислений с помощью программной платформы NVIDIA CUDA 117
4.4 Оптимизация программы вычисления с распараллеливанием по критерию машинного времени 120
4.4.1 Использование разделяемой памяти 121
4.4.2 Использование cuda аналогов математических функций 121
4.4.3 Выбор размера ячейки разбиения пространственной сетки области моделирования 122
4.5 Преобразования функции соударения сплайна и частницы инструментальной среды 124
4.6 Конфликты доступа процессов к глобальной памяти 127
4.7 Результаты сокращения времени моделирования 128
4.7.1 Сокращение времени в зависимости от количества частиц инструментальной среды и сегмента сплайна 128
4.7.2 Сокращение времени в зависимости от типа видеокарты 129
4.8 Моделирование стыкового профиля нервюры крыла ИЛ-96-300М 131
4.9 Оценка достоверности моделирования стыкового профиля нервюры крыла ИЛ-96-300М 134
4.10 Выводы по четвертой главе 141
Заключение 142
Список литературы 144
- Обзор методов моделирования многомассовых виброударных систем
- Описание основных этапов моделирования двумерных виброударных систем
- Алгоритм компьютерного моделирования процедуры интегрирования с распараллеливанием
- Оценка достоверности моделирования стыкового профиля нервюры крыла ИЛ-96-300М
Обзор методов моделирования многомассовых виброударных систем
Одним из первых методов моделирования динамики многомассовых систем в 1990-е годы был метод имитационного моделирования, базирующийся на их групповых свойствах [44, 68, 85].
При рассмотрении данного метода моделирования, интенсивно вибрирующие частицы (инструментальная среда) сначала рассматривалась как материальная точка [5], затем как упруго-вязко-пластичное тело [68], а в затем в виде послойной сплошной среды. Имитационное моделирование осуществляется в отсутствии информации о параметрах технологической системы, далее согласно итогам сопоставления теоретических характеристик, с экспериментальными, математическая модель наделяется упругодиссипативными и иными свойствами, обеспечивающими максимальное их совпадение [2]. Метод никак не предусматривает конфигурацию элемента и обладает в данном отношении неудовлетворительную достоверность. Точность изучения при имитационном моделировании по технологическим параметрам достигает более 100 %.
Моделирование методом интегральных оценок параметров инструментальной среды, представленной средой с накопленными (интегральными), упруго-диссипативными, зазорными и массовыми свойствами [44]. Осуществляется математическое моделирование основываясь на экспериментально определенных интегральных квазиупругих и диссипативных свойствах, динамических зазоров инструментальной среды, которые записывают в матрицу динамических свойств и методом перебора при постоянной частоте и фиксированных ускорениях, позволяет точнее определить динамические параметры процесса (фазовый угол соударения, импульс, силы и т.д.). Расчётные значения шероховатости, наклёпа и остаточных напряжений присваиваются всем участкам поверхности детали. Для деталей простой формы [43] погрешность метода составляет от 50 до 100 %. Существенным недостатком рассматриваемого метода является необходимость для каждой моделируемой детали опытного определения параметров инструментальной среды, что сопряжено с внушительными расходами. С целью уменьшения расходов в осуществлении экспериментальных исследований квазиупругие и диссипативные свойства инструментальной среды аналитически выражаются посредством динамических зазоров, которые проще определить экспериментально с помощью съемки цифровой камерой через прозрачную торцевую стенку контейнера [68]. Время моделирования составляет от 10 до 20 мин, но из-за учета лишь габаритных размеров детали и контейнера погрешности моделирования составляют пятьдесят – сто процентов.
В 2000-х годах, в связи с внезапным увеличением производительности персональных ЭВМ, возникли способы изучения, базирующиеся на прямом компьютерном моделировании с помощью событийно управляемых алгоритмов и методами молекулярной динамики [44, 75 – 81].
Использование метода моделирования конечного множества частиц с использованием метода дискретных элементов позволяет, кроме усреднённых значений динамических и технологических характеристик процесса, установить закономерности их распределения на различных областях элемента (детали) в поперечном сечении.
Определяются съём металла, шероховатость, остаточные напряжения и т.д. Выявляются координаты расположения недостаточно упрочнённых участков
Метод дискретного элемента (МДЭ) – это семейство численных методов специализированных с целью перемещения значительного числа элементов, базирующийся на приложении законов Ньютона и контактной механики, не имеющий недостатков непрерывных моделей, проявляющихся в результате дискретности его внутренней структуры [107]. Метод был поначалу использован Cundall в 1971 для решения вопросов исследования движения горных пород [93]. Затем в 1985 году Hocking, Mustoe, Williams совместно детализировали теоретические основы метода и показали, что МДЭ может рассматриваться, как совокупность метода конечных элементов (МКЭ, МДЭ). Применение данного метода к решению геомеханических задач описано в книге [100]. Williams, и Bicanic разместили несколько журнальных заметок описывающих современные направления в сфере МДЭ [105]. В книге [86] описано комбинирование МКЭ и МДЭ.
Суть МДЭ в том, что дискретный материал представляется совокупностью из N обособленных упругих частиц сферической формы диаметра D;. Движение i-той частицы, согласно законам классической Ньютоновской механики, полностью определяется двумя основными характеристиками: координатой центра тяжести рассматриваемой частицы Ц и углом её поворота вокруг центра тяжести фі.
Для каждой рассматриваемой частицы система уравнений её перемещения в прямоугольной декартовой системе координат может быть представлена в виде (1.1):
В свою очередь поверхностные силы Fy состоят из сил отталкивания и сил трения.
Методы дискретного элемента требуют большой объём вычислительных ресурсов ЭВМ, что ограничивает размер рассматриваемой модели или количество используемых частиц. Моделирование многомассовых виброударных систем МДЭ осуществляется в сечениях, совпадающих с плоскостью колебаний, при этом инструментальная среда представляется большим конечным множеством элементов (частиц) [44]. Геометрическая модель создается в масштабе, площадь сечений разбивается на маленькие области (ячейки), загружается элементами (частицами) инструментальной среды.
МДЭ применяемый для многомассовых виброударных систем позволяет определить значения удельной энергии соударений между частицами, энергии соударений границы упрочняемого объекта с частицами, установить съём металла, шероховатость, остаточные напряжения, наклёп на различных областях элемента (детали).
Расхождение теоретических и опытных результатов не превосходит значений от 40 до 50 %. Методы, разработанные ранее, имеют в два-три раза больше расхождения.
Моделирование виброударного упрочнения деталей с закреплением в контейнере осуществляется следующим образом [46-48, 61]. Формируется трёхмерная модель элемента (детали), контейнера и инструментальной среды; производятся поперечные разрезы в плоскости, совпадающей с траекторией колебаний, и продольные сечения; их контуры разбиваются на различные, прямолинейные отрезки. Площадь сечения делится на ячейки. В поперечных сечениях выполняется моделирование периодических перемещений указанных отрзков и центров частиц инструментальной среды, по которым определяются энергии периодических соударений между частицами и отрезками детали с частицами. Далее при помощи интерполяции строится пространственное распределение полученных технологических параметров [45].
В [33-42] создана модель трехмерной многомассовой системы и способ её изучения, обеспечивающие возможность параллельного расчёта перемещений элементов данной системы, а также их взаимодействий между собой. Автором [33] изобретен новый комплекс компьютерных программ для трёхмерного сеточного изучения многомассовых систем с отдельно распределёнными параметрами, осуществляющий человеко-машинный интерфейс, позволяющий изучить закономерности процесса с сохранением полученных данных в БД, а также быстрый просмотр и интерактивное представление цветовой информации в различных вариантах и разрезах без остановки процедуры моделирования.
Описание основных этапов моделирования двумерных виброударных систем
Охарактеризуем каждый из этапов моделирования. Отметим, что в той или иной мере все они были реализованы в работах [43-48, 56, 75-82]. Однако компьютерное моделирование на основе МДЭ приводит к необходимости выполнять большую вычислительную работу, требующую значительных затрат машинного времени (порядка 10 – 15 часов), что делает построенную модель технологически неэффективной. Основной задачей настоящей работы является оптимизация моделирования и организации вычислительного процесса с целью доведения их компьютерной реализации до приемлемых временных значений. Как показал предварительный анализ, проведенный автором настоящей диссертации и представленный в главе 3, наиболее трудозатратными по машинному времени являются этапы 3 – 6 (двумерного моделирования в плоскости), поэтому уделим больше внимания именно им.
Построение трехмерной геометрической модели системы контейнер-деталь инструментальная среда на первом этапе изображенного на рисунке 2.4, выполняется в трехмерном редакторе твердотельных конструкций, таком как Компас-3D – системы автоматизированного проектирования, разработанной российской компанией «АСКОН» с возможностями оформления конструкторской документации согласно стандартам серии ЕСКД или простого и производительного приложения для создания чертежей и 3D моделей любой сложности SOLIDWORKS 2017 [45,61].
В этих же прикладных пакетах выполняется второй этап построение набора характерных поперечных и продольных сечений модели системы, проиллюстрированных на рисунках 2.5, 2.6 и 2.7. В соответствии с полученными геометрическими параметрами построенных сечений – аналогичные сечения строятся (прорисовываются), как показано на рисунке 2.10 в программе моделирования [45,82].
На рисунке 2.11 дан поперечный схематический разрез вибросистемы
В этих же прикладных пакетах выполняется второй этап построение набора характерных поперечных и продольных сечений модели системы, проиллюстрированных на рисунках 2.5, 2.6 и 2.7. В соответствии с полученными геометрическими параметрами построенных сечений – аналогичные сечения строятся (прорисовываются), как показано на рисунке 2.10 в программе моделирования [45,82]. деталью и контейнером, заполненным инструментальной средой. В дальнейших рассмотрениях расчётная область, изображённая на рисунке 2.8 – круг (прямоугольник) с лакуной (сечением детали) внутри, который заполняется мелкими абразивными элементами (круглыми частицами). Внешние воздействия моделируются функциями, задающими гармонические колебания границы (контуры контейнера и детали) вдоль координатных осей.
Затем происходит разбиение аппроксимированной границы на отрезки сплайнов, представляющие собой кусочки (прямых, парабол, эллипсов) примерно одинаковой длины и осуществляется вычисление координат этих отрезков в выбранной и фиксированной системе координат. За длину сплайна принимается радиус самой мелкой частицы в ансамбле частиц
Направление отрезков на плоскости устанавливает ориентацию нормали создаваемых сплайнов, которая в том числе показывает направление размещения элементов системы. Результатом вычислительной работы на этом этапе является определение и фиксация начального положения и ориентации сплайнов на плоскости (от 2000 до 3000 сплайнов).
Загрузка контейнера (рисунок 2.8) инструментальной средой, осуществляемая на четвертом этапе, производится в следующем порядке. Вычисляется количество частиц инструментальной среды согласно формуле
Моделирование и численный расчёт временных характеристик виброударного упрочнения (в плоскости) на основе метода дискретных элементов(МДЭ) реализуется на пятом этапе. При этом инструментальная среда представляется набором из N круглых упругих частиц диаметра D,. Обозначим ft - радиус-вектор /-ой частицы, йі - её скорость, mt - массу, ф1- угол поворота вокруг центра тяжести, /. - момент инерции и сдг- угловую скорость. Движение (поступательное и вращательное) каждой частицы полностью описывается основными уравнениями динамики
Расчёт контактных сил может проводится с учётом различных теоретических моделей вязкоупругого взаимодействия частиц. В настоящей работе используется вязкоупругая модель Кельвина-Фохта, успешно работающая в случаях, когда статическая характеристика нагружения в контакте близка к линейной. Учет силы трения является одним из наиболее важных вопросов в рассматриваемой задаче. В разработанной программе моделирования реализованы следующие модели взаимодействия между частицами инструментальной среды и частицами с контейнером и деталью: кулоновское трение, использующее классический закон Кулона, контакт Уолта-Людинга[56], модели вязкого трения Фойхта-Кельвина, модели Финни-Рабиновича [56]. Подробное описание построения контактных сил можно найти в [56], [82]. Мы будем пользоваться полученными в [56] формулами. Отметим, что с технологической точки зрения необходимо рассматривать два различных случая:
1) взаимодействие частиц инструментальной среды между собой;
2) взаимодействие отрезков сплайнов контейнера и детали с частицами инструментальной среды.
Алгоритм компьютерного моделирования процедуры интегрирования с распараллеливанием
Рассматривая алгоритм последовательного компьютерного моделирования процедуры интегрирования двумерной многомассовой системы методом Адамса 4 порядка с предикт-коррекцией (шаги 14-27, рисунок 2.22) выявлено, что, во-первых, самыми длительными по времени являются шаги (18-21, 24-26), так как именно здесь происходит расчет сил, действующих на все элементы системы и между ними и, во-вторых, вычисления на шагах 14-17, 22-23 и 27 можно легко распараллелить. Так как количество элементов системы лежит в диапазоне от нескольких сотен до нескольких десятков тысяч, то необходимо при распараллеливании вычислений применить аппаратные средства, которые могут генерировать примерно такое же количество потоков. В качестве таких средств в настоящей работе используется видеокарты с поддержкой технологии NVidia CUDA.
Процедура интегрирования двумерной многомассовой системы с распараллеливанием выполняется на основной системе (называемый термином хост), к которой подключен GPU – сопроцессор для массивно-параллельных вычислений, который по отношению к хосту называют устройством. CUDA – программа задействует CPU, где выполняется последовательная часть кода и подготовительные стадии для GPU вычислений, а параллельные участки кода переносятся на GPU, где будут выполняться большим множеством нитей.
Исходный код программы прямого моделирования реализован в среде Delphi. Для того, чтобы применить распараллеливание при помощи NVIDIA CUDA необходимо текст кода с распараллеливанием операций реализовать в среде C++. Cвязано это с тем, что на текущий момент времени компилятор CUDA поддерживает только языки С, С++ и Fortran для ускорения приложений с помощью графических процессоров NVIDIA с массивно параллельной архитектурой. Все функции, написанные на С++ реализовываются в форме динамической библиотеки, которая подключается к исходной программе моделирования на языке Pascal.
Алгоритм моделирования с распараллеливанием продемонстрированный на рисунке 3.11 отличается тем, что процедура инициализации, которая производится вызовом функций из динамической библиотеки, выполняется как в общей ОЗУ хоста, так и устройства. Затем после выполнения процедуры разгонки «warmup» шаг 13 рисунка 2.23, выполняющейся на CPU осуществляется копирование всех данных, необходимых для моделирования системы (координат всех отрезков сплайнов, координаты и скорости частиц инструментальной среды и т.д.) с хоста на устройство и последующие 100 итераций интегрирования выполняются на устройстве посредством вызова функций из динамической библиотеки. Каждая из вызываемых функций реализована в виде вызова ядер CUDA, код в которых выполняется параллельно, т.е. внутри ядра одновременно происходят вычисления для каждой из частиц инструментальной среды и каждого сегмента детали. Таким образом, распараллелены вычисления на наиболее длительных шагах 14-27 показанных на рисунке 2.23.
Алгоритм моделирования с распараллеливанием включает следующие этапы:
Шаг 1. Инициализация всех данных в общей ОЗУ.
Шаг 2. Инициализация устройства GPU включающая в себя:
выбор наибольшего по производительности устройства GPU,
резервирование памяти GPU под структуры данных, которые в последующем будут использоваться в программе.
Все указанные действия производятся вызовом функций из динамической библиотеки.
Шаг 3. Процедура разгонки «warmup» – выполняется три шага интегрирования системы явным двухшаговым методом. Данная процедура выполняется на CPU.
Шаг 4. Копирование всех данных, необходимых для моделирования системы (координат всех отрезков сплайнов, координаты и скорости частиц инструментальной среды и т.д.) с хоста на устройство.
Шаг 5. Обнуление счетчика итераций.
Шаг 6. Если счетчик итераций не достиг максимального количества шагов, то увеличиваем его на 1, если достиг, то переходим к шагу 13.
Шаг 7. Выполняется 1 шаг интегрирования на ГП посредством функций из динамической библиотеки.
Шаг 8. Если на данном шаге требуется сохранить данные, то копируем данные из ГП на хост, если нет, то переходим к шагу 10.
Шаг 9. Процедура сохранения данных.
Шаг 10. Если на данном шаге требуется визуализация, то копируем данные из ГП в хост посредством функций из динамической библиотеки, если нет, переходим к шагу 6.
Шаг 11. Визуализация текущих значений фазовых траекторий, которая реализована функционалом класса TCanvas визуальных компонентов Delphi.
Шаг 12. Переход к шагу 6.
Шаг 13. Копирование данных из памяти ГП в общую память. Шаг 14. Сохранение данных на диск.
Шаг 15. Освобождение выделенной памяти под структуры данных на ГП.
Шаг 16. Освобождение выделенной памяти под данные на хосте.
Шаг интегрирования с распараллеливанием полностью выполняется на GPU и состоит из функции ядра, которые выполняются по одной независимой нити для каждого набора элементов и вспомогательных device функций.
Алгоритм одного шага интегрирования с распараллеливанием вычислений изображённый на рисунке 3.12 состоит из следующих этапов:
Шаг 1. На первоначальном этапе процедура интегрирования принимает значение времени time = time + step, где step - шаг по времени.
Шаг 2. Выполнение ядра copyBoundary - копирование всех параметров сплайновой границы из предыдущего состояния в текущее.
Шаг 3. Получаем состояние сплайновой границы в момент времени time , вызывая ядро buildboundary . На данном этапе происходит копирование параметров отрезков сегмента сплайна (координат начала и конца, длины, нормали и т.д.), параметров пчевдочастицы (позиции, скорости, угловые скорости, фазы), технологических параметров отрезка сплайна (начала контакта, продолжительность контакта, скорость соударения, съем, наклеп, остаточные напряжения и т.д.).
В завершении текущего этапа, в зависимости от выбранного закона движения сплайновой границы меняем координаты и скорости отрезка сплайна и псевдочастицы.
Шаг 4. Ядро predictor - ядро, реализующее схему предикт-коррекции для положений, скоростей, угловых скоростей и фаз ансамбля частиц инструментальной среды.
Шаг 5. makePosVel - вспомогательная функция ядро, осуществляющее копирование позиций и скоростей частиц из двумерного массива в одномерный. Одномерный массив далее используется в процедуре индексирования.
Шаг 6. Ядро calcHash вычисляет хэш-значение для каждой частицы на основе ее идентификатора.
Шаг 7. Ядро sortParticles осуществляет сортировку частиц на основе их хэш-значения. Сортировка производится с помощью библиотеки thrust – одной из самых красивых и гибких библиотек для CUDA. Основным отличием данной библиотеки от других распространенных библиотек для CUDA является то, что thrust - это библиотека, основанная на использовании шаблонов (template) языка С++.
Все классы и функции в этой библиотеки – шаблонные, т.е. для использования данной библиотеке необходимо только лишь подключить соответствующие заголовочные файлы. Thrust содержит набор различных параллельных примитивов, таких как различные преобразования, сортировка, операции reduce и scan. Все последние версии CUDA включают в себя thrust, а это значит, что для работы с thrust никаких дополнительных установок не требуется.
Шаг 8. Ядро findCellStart поиск начала любой заданной ячейки в отсортированном списке. Текущий код также находит индекс в конце каждой ячейки аналогичным образом.
Шаг 9. Функция D, выполняемая на GPU содержит три ядра:
Ядро D_Device – расчёт взаимодействий между частицами инструментальной среды,
Ядро D_Spl_Device – расчёт взаимодействий между частицами инструментальной среды и сплайновой границей,
Ядро Overrading – осуществляет пересчёт позиций и скоростей и фаз с учетом полученных сил в предыдущих ядрах.
Шаг 10. Ядро corrector реализует схему корректора для положений, скоростей, угловых скоростей и фаз ансамбля частиц инструментальной среды.
Шаг 11. Вызов функции D.
Шаг 12. Ядро swapArrays, выполняющее процедуру присваивания – шаг 27 рисунка 2.22.
Шаг 13. Переходим к шагу 1.
Оценка достоверности моделирования стыкового профиля нервюры крыла ИЛ-96-300М
Моделирование стыкового профиля нервюры крыла ИЛ-96-300М в последовательной программе с экспериментальной проверкой достоверности результатов моделирования подробно проводилось в работе [61] с участием Копылова Ю.Р и автора данной диссертации.
В данном разделе представлены графические результаты моделирования процесса виброударного упрочнения стыкового профиля нервюры крыла ИЛ-96-300М рисунки 4.14 – 4.19 в последовательной реализации программы в среде Delphi [61] в сравнении с параллельной реализацией ввиде комплекса программ, использующего NVIDIA CUDA, описание которого содержится в данной диссертации.
Траектория движения центра масс конечного множества частиц в последовательной а) и параллельной б) программе представлена на рисунке 4.14.
Длина пробега центра масс конечного множества частиц, характеризующая амплитуду динамического зазора между частицами, в а) последовательной и б) параллельной программе представлена на рисунке 4.15.
Графическое изображение модуля скорости виброударного перемещения частиц инструментальной среды в последовательной (а) и параллельной программе (б) представлены на рисунке 4.17.
Гистограмма диапазонов скоростей частиц инструментальной среды в последовательной (а) и параллельной программе (б) представлены на рисунке 4.18.
Обоснован выбор дизайна программы моделирования, использующий взаимодействие программных платформ Delphi и C++ в рамках распараллеливания вычислительных задач.
Разработан новый комплекс программ для моделирования двумерных многомассовых систем с распараллеливанием вычислений при помощи программной платформы NVIDIA CUDA.
Описана оптимизация программы двумерного моделирования многомассовых систем.
Представлены результаты сокращения времени моделирования в зависимости от количества частиц инструментальной среды, количества сегментов сплайна, а также типа видеокарты.
Представлены результаты моделирования стыкового профиля нервюры крыла ИЛ-96-300М.
Представлены результаты сокращения времени моделирования стыкового профиля нервюры крыла ИЛ-96-300М: при моделировании на видеокарте Nvidia GeForce 970 достигается 25-ти кратное ускорение вычислений (26 мин) в сравнении с последовательной реализацией программы, которая выполняется в течение 660 мин.