Параллельный алгоритм реконструкции поверхности прочности композиционных материалов для архитектуры Intel Mic (Intel Many Integrated Core architecture) Текст научной статьи по специальности «Физика»

УДК 004.272.2

А. П. Соколов, В. Н. Щетинин, А. С. Сапелкин

Параллельный алгоритм реконструкции поверхности прочности композиционных материалов для архитектуры Intel MIC (Intel Many Integrated Core Architecture)

Аннотация. Целью исследования было создание параллельной программной реализации численного метода реконструкции моделей поверхности прочности первичного разрушения исследуемых композиционных материалов. Использовался квадратичный критерий прочности Малмейстера—Ву. В основе использовались методы асимптотического осреднения (Бахвалов Н. С., Победря Б. Е.) и конечных элементов. Программная реализация была создана в рамках графоориентированной технологии, реализованной для архитектуры Intel MIC в Распределенной вычислительной системе GCD. Были проведены вычислительные эксперименты для серии моделей композиционных материалов, задаваемых их схемами армирования («ячейками периодичности»), по определению сечений поверхности прочности первичного разрушения. Рассматривались SD-армированный и ID-армированный композиты. Результаты расчетов представлены.

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

Введение

В процессе решения прикладных задач анализа механических характеристик композиционных материалов (КМ) часто применяют численные подходы, одним из которых является метод асимптотического осреднения (МАО), называемый в зарубежной литературе гомогенизацией [1,2]. В основе метода лежит решение серии краевых задач Ьрд специального вида, называемых «локальными» [3], каждая из которых решается методом конечных элементов (МКЭ) [4,5].

Работа проведена на инициативной основе. © А. П. Соколов, В. Н. Щетинин, А. С. Сапелкин, 2016 (с МГТУ имени Н.Э. Баумана, 2016 © Программные системы: теория и приложения, 2016

МАО обеспечивает возможность расчета эффективных характеристик (ЭК) периодических гетерогенных структур, с помощью которых могут быть описаны схемы армирования композиционных материалов, используемых в различных конструкциях. К таким характеристикам следует отнести упругие (тензоры модулей упругости cijki), тепловые (тензоры теплопроводности «¿¿), электростатические, прочностные (пределы упругости, текучести) и прочие характеристики. После определения ЭХ композитные элементы исследуемой конструкции могут быть заменены на эквивалентные однородные с полученными свойствами. Далее, к примеру, в процессе конечно-элементного решения задач о напряженно-деформируемом состоянии (НДС) композитные элементы конструкции могут быть разбиты на значительно более грубую вычислительную сетку, что позволяет существенно снизить вычислительные затраты.

Объектом настоящего исследования были эффективные упруго-прочностные характеристики (ЭУПХ) композиционных материалов, описываемых различными схемами армирования.

Задача расчета ЭУПХ КМ при фиксированной нагрузке с использованием МАО требует существенных вычислительных ресурсов.

Для определения сечения поверхности прочности могут потребоваться многократные вычисления ЭУПХ КМ для различных путей нагружения в плоскости, заданной двумя выбранными компонентами тензора напряжений ст^. Такой расчет для нескольких сотен путей нагружения приведет к часовым временным затратам.

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

В настоящей работе был разработан параллельный алгоритм вычисления пределов первичного разрушения ар^mnpst . Программная реализация была построена с использованием графо-ориентированной технологии [6], применяемой для создания программных реализаций сложных вычислительных методов в рамках Распределенной вычислительной системы GCD. Вычислительные эксперименты были проведены с использованием высокопроизводительной вычислительной системы на основе архитектуры Intel MIC, что позволило получить сечения поверхностей прочности серии различных композиционных материалов.

1. Постановка задачи построения поверхности прочности композиционного материала

1.1. Прочностные характеристики композиционного материала для пространственного случая

Напряженно-деформированное состояние материала в общем случае определяется тензором напряжений а. Для пространственного случая в классической теории упругости тензор напряжений симметричен и имеет 6 независимых компонент. Существует ряд математических теорий прочности [7-12], и все они базируются на введении и вычислении специальной скалярной функции повреждаемости с последующим ее сравнением с единицей:

(1) F(aij(x),Sa(x)) V 1,

где x — декартовы координаты точки,

aij (x) - текущее напряженно-деформированное состояние в точке x, характеризуемое значениями компонент тензора напряжений, Sa(x) - прочностные характеристики материала в точке x.

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

(1) ввести универсальное обобщенное обозначение технических пределов прочности, а именно пределов пропорциональности, текучести и полного разрушения материала в точке;

(2) выбрать наиболее универсальный критерий прочности материала, который мог бы быть применим для большинства материалов.

Обозначение 1.1. Введем общее обозначение технических пределов прочности (пределов пропорциональности, текучести и полного разрушения) в 6-мерном пространстве тензора напряжений:

atLMNPST,

где t € {р, у, и};

где р - обозначение для пределов упругости или пропорциональности (Proportional limit);

у - обозначение для пределов текучести (Yield strength); и - обозначение для пределов полного разрушения (Ultimate strength); {L, М, N, Р, S, Т} € {0, Е, С}, где Е - индикатор растяжения (Extension), С - индикатор сжатия (Compression);

0 - указывает на отсутствие НДС по соответствующей компоненте тензора напряжений;

положение каждого индекса соответствует определенной компоненте тензора напряжений: L ^ оц, М ^ о"22, N ^ <733, Р ^ (Jyi, $ ^ ^13, Т ^ ^23-

Например, используя введенные обозначения: аеЕ00000 = — предел пропорциональности при растяжении вдоль оси ОХ; аиоосооо = ас — предел полного разрушения при сжатии по оси OZ; аи000Е00 = as — предел полного разрушения при сдвиге в плоскости OXY; ауЕ00000 = at — предел текучести при растяжении вдоль оси ОХ.

Обозначение 1.2. Простым путем нагружения будем называть прямую, индуцированную базисным вектором = Pi~£i =

(pi,P2,P3,P4,P5,Pe)T, i € {1, 2, 3, 4, 5, 6} определенным в 6-мерном пространстве тензора напряжений, где pi ^ ац, Р2 ^ 022, Р3 ^ &33, Ра ^ &12, Ръ ^ ^13, Рв ^ °23;

Pi € {0,1, —1}, i € {1, 2, 3, 4, 5, 6}, где значение '1' его компоненты указывает на увеличение (или растяжение), а значение '-1' его компоненты указывает на уменьшение (или сжатие) соответствующей компоненты тензора напряжений.

Обозначение 1.3. Будем называть эффективными техническими упруго-прочностными характеристиками исследуемого КМ:

(1) технические упругие постоянные, получаемые на основе компонент эффективного тензора модулей упругости Cijki и эффективного тензора упругих податливостей П-ijki при допущении об ортотропности модели гомогенизированного композиционного материала:

• модули Юнга Ex,Еу,Ez;

• коэффициенты Пуассона vxy , vyz , vxz;

• модули сдвига Gxy, Gyz, Gxz;

(2) осредненные технические пределы прочности первичного разрушения в обозначениях, введенных выше:

• lmnpst, где {L,M,N,P,S,T} € {0,Е,С}.

• параметры L, М, N, Р, S, Т определяются выбранным в данный момент простым путем нагружения (см. Обозначение 1.2).

Следует отметить, что общее число М простых путей нагружения, исходя из введенного Обозначения 1.2, определяется по формуле

М = 36 — 1 = 728 [7], что фактически определяет общее число пределов прочности, которые могут быть найдены для выбранного КМ в рамках данной теории прочности.

1.2. Упрощенный квадратичный критерий прочности

Будем использовать критерий прочности, записанный в следующей форме:

(2) Ва + 3 Вха2 + В2а2 = 1.

Константы В, В\, В2 можно выразить через пределы прочности ат на растяжение, ас — на сжатие, а8 — на сдвиг:

(3) В = 1 — ±,

ат ас

(4) В1 3 1

атас a2s

(5) в2 = 3-L,

где а = 1\((Тц) - первый инвариант тензора напряжений, а^ - интенсивность тензора напряжений:

(6) а, = sj ((on - ^22)2 + (0-11 - озз)2 + {oil - озз)2 + 6(gf2 + of3 + о-2з))

для декартовой системы координат,

(7) h(a) = an дгг = ап + а2 2 + о33. 1.2.1. Инверсия прочности

Случай, когда ас = ат ас = av, т.е. когда прочность материала на сжатие и растяжение одинакова, называют инверсией прочности [7]. В этом случае коэффициент В = 0.

Не ограничивая общности в формулах при этом используют только предел прочности на растяжение ат. С учетом этого факта преобразуем (2), подставив в соотношение критерия прочности выражения для констант прочности В\,В2, получим:

К ^ - £) + ¿2 '2 = 1

К1 - ^22)2 + К1 - ^33)2 + (^22 - ^33)2 + + а13 + а2з)

6а2

+

+ ( "2 - Б"? ) (а11 + а22 + ^33)2 = 1.

Раскроем скобки и упростим полученное соотношение. Получаем вид упрощенного квадратичного критерия прочности:

(8)

2 2 2 2 (о~11 + ^22 + ^зз) + ^12 + а13 + ^3 - ^11^22 - ^11 ^33 - ^22^33 = -

2 + 2 1.

1.2.2. Устранение допущения об инверсии прочности

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

Для этого следует различать прочностные характеристики на сжатие ас и растяжение ат. Тогда представленный выше критерий прочности принимает вид:

(9) - ^г) (ои + (?22 + (У33) + 4т(^11 + а22 + ^33)2 +

1 ( 2 2 2 \ 1

+ 412 \а12 + а13 + аЖ - а11а22 - ^11^33 - ^22^33) = 1.

1.3. Квадратичный критерий прочности в общем виде

Наиболее широко используемым критерием прочности является тензорно-полиномиальный критерий Гольденблата-Копнова или в некотором частном случае Цая-Малмейстера [7].

Введем функцию повреждаемости, как функцию компонент тензоров прочности и напряжения, следующим образом:

(10) ^ (в(1),3(2),а) = + в™ , аг,ак1,

где Я(1), Я(2) - тензоры прочности второго и четвертого рангов соответственно, определяемые прочностными свойствами материала, а а -тензор напряжений.

Критерий прочности записывается так:

(11) ^ (в (1),в(2) , а) = 1.

При выполнении этого равенства достигается разрушение.

1.4. Реконструкция поверхности прочности

Обозначение 1.4. Поверхностью прочности называется геометрическое место точек заданных неявно уравнением (1), с явно заданной функцией повреждаемости F(aij(x), Sa(x)).

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

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

1.5. Формальная постановка задачи

Исходные данные.

(1) Модель композиционного материала: задана характерным объемом области КМ, определяющим схему армирования всего материала, и называемым ячейкой периодичности (ЯП) V. ЯП КМ V = У Va, где a G N, а — номер компоненты КМ (обычно

а

выделяют наполнитель а = 1 и матрицу а = 2) [13]. Ячейка может обладать многоуровневой структурой, — в таком случае задано древовидное описание микроструктуры КМ [14], для которой определены геометрии всех необходимых ЯП для каждого структурного уровня где п — номер уровня. Геометрия ЯП

задавалась с применением технологии Constructive Solid Geometry (CSG). Каждый узел древовидного (или иерархического) описания микроструктуры КМ должен быть представлен своей ячейкой периодичности и для каждой ЯП должен быть определен набор материалов-компонентов.

(2) Для каждого материала-компонента каждой ячейки периодичности иерархического описания микроструктуры КМ должны быть заданы упруго-прочностные свойства:

(a) минимально достаточный набор упругих констант для определения компонент тензора модулей упругости С<а^ы: модули Юнга Ex ,Еу ,Ez ; коэффициенты Пуассона vxy , vyz, vxz; модули сдвига Gxy , Gyz , Gxz;

(b) минимально достаточный набор пределов прочности (пропорциональности) apLmnpst каждого материала-компонента КМ для определения компонент тензоров прочности.

Необходимо определить:

(1) компоненты эффективных тензоров модулей упругости Сцк1 и упругих податливостей П-до исследуемого КМ;

(2) эффективные упругие технические константы КМ: модули Юнга, коэффициенты Пуассона, модули сдвига;

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

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

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

Метод асимптотического осреднения был предложен проф. Бахва-ловым Н.С. и проф. Победрей Б.Е. [13,15] как инструмент, обеспечивающий математическое описание гетерогенных сред с периодической микроструктурой.

Подробно теория МАО и его программная реализация рассмотрены в диссертационном исследовании [16].

В рамках настоящего исследования МАО использовался как основа для расчета прочностных характеристик исследуемых КМ.

1.7. Алгоритм поиска предела пропорциональности КМ для выбранного простого пути нагружения

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

Аналитический метод определения предела пропорциональности.

(1) Зададим простой путь нагружения " = ар^ "г, где а € И — длина вектора, "" — базисные векторы в декартовой системе координат в 6-ти мерном пространстве тензора напряжений, г €{1,..., 6}. Компоненты р^ этого базисного вектора не могут задаваться независимо, так как для прикладных задач всегда нужно определять пределы прочности по определенным осям или в определенных плоскостях. Поэтому в большинстве случаев имеем, к примеру: (р1,Р2,Р3,Р4,Р5,Ре)Т = (1,0, 0, 0,0, 0)Т

— базисный вектор простого пути нагружения при растяжении по оси ОХ, (р1,р2,рз,р4,рв,р6)т = (0,0, —1, 0,0, 0)Т — базисный вектор простого пути нагружения при сжатии по оси OZ, (р\,р2,рз,Р4,Р5,Рб)Т = (0,0, 0,1, 0,0)Т — базисный вектор простого пути нагружения при чистом сдвиге в плоскости ОХУ и т.д. Таким образом компоненты становятся зависимыми и при увеличении или уменьшении нагрузки меняются синхронно.

(2) Тензор концентрации напряжений связывает средние напряжения аыI с локальными микронапряжениями в ЯП, что верно только для линейно-упругой модели деформирования. Используя ¿ы (£) могут быть найдены поля компонент тензора микронапряжений а^(^) = в ячейке периодичности без полного решения разрешающей СЛАУ МКЭ.

(3) Простой путь нагружения , «домноженный» на действительный параметр а, определяет некоторые средние напряжения а^. При некотором а* достигается первичное разрушение материала.

(4) Подставив средние напряжения и пределы прочности конкретной компоненты в уравнение

(12) ^(аВцЫ(€)ам,в, Я(2)) = 1,

будет получено уравнение относительно а* (для квадратичного критерия Малмейстера-Ву квадратное уравнение). Для несложных критериев прочности это уравнение может быть решено аналитически.

(5) Получив а* корень уравнения (12), может быть вычислена длина вектора , которая и будет искомым пределом прочности первичного разрушения.

2. Программная реализация вычислительного алгоритма

Разработка программной реализации алгоритма решения поставленной задачи была осуществлена с использованием языка программирования С+—Н и графоориентированной технологии построения программных реализаций сложных вычислительных методов [6] Распределенной вычислительной системы ССБ.

2.1. Особенности последовательной программной реализации алгоритма

(1) Инициализация объектов:

(а) загрузка конечно-элементной сетки (Mesh);

(b) выбор типа конечного элемента (Finite Element Type);

(c) определение фиксированного набора из 6 граничных условий (ГУ) для задачи поиска эффективных упругих характеристик, которые необходимы для применения метода гомогенизации (BoundaryConditionsMap) к каждой соответствующей задаче Lp q ;

(d) загрузка данных о материалах-компонентах исследуемого КМ (SLD — Solid data).

(2) Для всех краевых условий (для каждой задачи Lpq ) должна быть решена задача о напряженно-деформированном состоянии согласно следующей последовательности операций:

(a) инициализация объекта решателя разрешающей системы линейных алгебраических уравнений (СЛАУ) (LASSolver);

(b) ассемблирование матрицы жесткости и вектора правой части;

(c) решение СЛАУ;

(d) по найденному полю вектора перемещений Ui решению в узлах поиск компонент тензора деформации £ij и тензора напряжений ai j.

(3) Поиск компонент тензора концентрации напряжений Вц^.

(4) Расчет эффективных упругих свойств КМ.

(5) Расчет пределов первичного разрушения.

(a) Инициализация набора простых путей нагружения.

(b) Вычисление пределов первичного разрушения (пропорциональности) ap L MNPST.

2.2. Особенности параллельной программной реализации алгоритма

2.2.1. Технические характеристики вычислительного узла

Программа алгоритма вычисления точек сечений поверхностей прочности исследуемых КМ была реализована на вычислительном узле с архитектурой Intel MIC. Узел включал в себя: управляющий хост на базе 2-х процессоров Intel Xeon E5-2695 с общей памятью и сопроцессор Intel Xeon Phi 7120, подключенный к PCI слоту хоста. Технические характеристики вычислительного узла:

• Intel Xeon E5-2695: 12 физических ядер с тактовой частотой 2.4 GHz (24 логических ядра), объем оперативной памяти 64 Гб, объем кэш памяти ядра 30 Мб, длина векторного регистра 256 бит.

• Intel Xeon Phi 7120: 61 физическое ядро с частотой 1.238 GHz каждое (244 логических ядра), объем оперативной памяти 16 Гб, объем кэш памяти ядра: 30 Мб, длина векторного регистра 512 бит.

Для компиляции был использован специализированный программный пакет Intel Parallel Studio 2016 XE Composer Edition.

2.2.2. Параллельный алгоритм поиска тензоров концентраций напряжений и упругих податливостей КМ

Легко заметить, что основная вычислительная часть последовательной схемы алгоритма может быть легко распараллелена: задачи Lpq о НДС для разных ГУ независимы и могут быть решены параллельно. Выделим пункты 1-3 последовательной схемы (см. стр. 11) в отдельный алгоритм поиска тензоров концентрации напряжений Bijki и упругих податливостей П^ы КМ для возможности исследования эффективности его распараллеливания на процессоре Intel Xeon E5-2695. Суть параллельной схемы решения задачи заключается в независимом решении задач Lpq с выделением для каждой отдельного вычислительного потока. Многопоточное программирование осуществлялось при помощи технологии OpenMP. Отметим, что алгоритм не потребовал специальной подготовки данных для каждого потока: все данные инициализированные в пункте 1 схемы являются разделяемыми.

Схема параллельного алгоритма поиска тензоров концентраций напряжений и упругих податливостей КМ:

(1) Инициализация объектов (аналогично последовательной схеме).

(2) Инициализация 6-ти потоков в соответствии с числом задач Lpq 1, где {p,q} G {1, 2, 3}.

(3) Решение задач Lpq МКЭ в отдельных вычислительных потоках.

(4) Барьерная синхронизация потоков по окончанию вычислений, слияние результатов вычислений.

(5) Расчет тензора концентрации напряжений B^ki (аналогично последовательному алгоритму).

(6) Расчет эффективных упругих свойств КМ (аналогично последовательному алгоритму).

1Для задачи поиска эффективных упругих характеристик КМ количество задач Ьрц шесть, тогда как для поиска эффективных тепловых характеристик количество задач будет 3. Таким образом количество требуемых потоков для наиболее эффективного расчета зависит от типа задачи.

Рис. 1. Сравнительные временные затраты при решении задачи на различных вычислительных сетках

(7) Расчет пределов первичного разрушения (аналогично последовательному алгоритму).

3. Вычислительный эксперимент

3.1. Сравнение ускорения при расчете на различных вычислительных сетках

Вычислительные эксперименты проводились на различных вычислительных системах, для которых далее проводился сравнительный анализ затрат времени на решение поставленной задачи методом гомогенизации. Обеспечивались: параллельное решение задач Lpq и параллельное решение систем линейных алгебраических уравнений (СЛАУ) в рамках одной Lpq методом сопряженного градиента (CG) c учетом использования оптимизирующих компиляторов Intel Compiler и дополнительных команд ассемблера AVX.

Используемые вычислительные системы:

• эталонная рабочая станция (РС);

• вычислительная серверная система на основе архитектуры Intel

MIC;

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

Таблица 1. Значения времени решения в секундах

Число узлов в КЭ сетке 12355 36150 300482

Последовательный алгоритм 20 98 1531

Параллельная обработка Lpq, Intel MIC 7.7 28.7 377

Параллельное решение СЛАУ методом CG, Intel MIC 7.0 21.6 277

После векторизации вычислительных операций, Intel MIC 2.4 10.7 140

3.2. Параллельное решение СЛАУ в рамках метода сопряженного градиента CG

Для решения разрешающей СЛАУ МКЭ использовался метод сопряженного градиента (Conjugate Gradient (CG)). Основная ресурсоемкая операция этого итерационного метода — вычисление произведения матрицы на вектор. Эта операция хорошо распараллеливается как на системах с общей памятью, так и на системах с распределенной памятью. В предыдущем пункте для решения задачи были задействованы только 6 вычислительных ядер процессоров хоста из 48 доступных.

Для распараллеливания метода CG были задействованы оставшиеся 42 логических ядра. Для хранения СЛАУ использовался формат CSR (Compressed Sparsed Row).

Заметим, что такой подход обеспечивает существенное ускорение только на задачах с большим числом узлов и малое — на задачах с малым числом узлов.

3.3. Использование векторных инструкций

Используемый процессор Intel Xeon E5-2695 имеет векторный регистр длиной 256 бит, что позволяет, используя набор команд ассемблера AVX, векторизовать следующие части алгоритма:

(1) решение СЛАУ;

(2) вычисление компонент тензоров деформации и напряжений;

(3) вычисление компонент тензора концентрации напряжений.

После векторизации были получены результаты, представленные на рис. 1 и в сводной таблице 1.

Прирост производительности при использовании команд ассемблера AVX в совокупности с параллельным решением задач Lpq и параллельным решением СЛАУ очевиден.

3.4. Параллельный алгоритм поиска пределов

пропорциональности по заданному набору простых путей нагружения

С точки зрения параллельных вычислений расчет отдельного предела прочности по заданному фиксированному простому пути нагружения никак не связан с аналогичным расчетом для другого пути, что позволяет эффективно задействовать архитектуру Intel MIC в части многопоточного распараллеливания вычислений по отдельным путям нагружения2.

3.4.1. Параллельный поиск пределов пропорциональности на хосте

Был проведен поиск пределов пропорциональности только на хосте. Так как общее число логических ядер равно 48, то общее число простых путей нагружения были поделены на группы по 48. Расчет по каждой группе был проведен в отдельном OpenMP потоке.

Расчеты проводились на тех же конечно-элементных аппроксимациях ячеек периодичности исследуемых композитов. Заметим, что одна из размерностей тензора В совпадает с количеством конечных элементов, соответственно по ней можно судить о вычислительной сложности задачи. Количество точек поверхности прочности для всех экспериментов было задано равным 10 000.

Результаты представлены на рис. 2, сравнение абсолютных значений затрат времени представлен в таблице 2).

Отметим, что с ростом вычислительной нагрузки в параллельной части алгоритма производительность также возрастает.

3.4.2. Использование технологии Offload передачи данных между хостом и сопроцессорами

Применение технологии Offload позволяет перенести вычислительную нагрузку с основного процессора Intel Xeon E5-2695 на сопроцессор Intel Xeon Phi 7120.

Результаты представлены на рис. 2, сравнение абсолютных значений затрат времени представлен в таблице 2).

2В ближайшее время планируется перенести эту вычислительную нагрузку на сопроцессор.

Рис. 2. Сравнение временных затрат при решении задачи реконструкции сечений поверхностей прочности исследуемых КМ при выборе различных путей нагружения.

Таблица 2. Значения времени в секундах

Число узлов в КЭ сетке 12355 36150 300482

Последовательный алгоритм, PC 850 4550 37041

Параллельная обработка

параллельная обработка отдельных 29 7729 102 802 1714 51 путей нагружения в отдельных потоках (Intel Xeon E5-2695)

Параллельная обработка отдельных

путей нагружения в отдельных 11.233 33.87 333.435

потоках на сопроцессоре Intel Xeon Phi 7120

3.5. Примеры графических результатов расчетов

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

(1) Ш-дисперсно-армированного композиционного материала;

(2) 3В-ортогонально-армированного композиционного материала.

(a) при сдвиге (b) при сжатии

Рис. 3. Решения задач Lpq в рамках метода гомогенизации: поля функции повреждаемости при деформации вдоль оси армирования

3.5.1. Исследование lD-дисперсно-армированного композиционного материала

Отдельные результаты решения промежуточных задач Lpq в рамках метода гомогенизации представлены на рис. 3. Необходимость расчета упруго-прочностных характеристик представленных материалов возникает, например, при проектировании специализированных газоразделительных модулей мембранного типа.

Полученные отдельные сечения поверхности прочности первичного разрушения lD-дисперсно-армированного композиционного материала в выбранных осях представлены на рис. 4(a), 4(b), 4(c).

3.5.2. Исследование Эй-ортогонально-армированного композиционного материала

3D-ортогонально-армированного композиционного материала (3D КМ) активно применяются при проектировании теплозащитных конструкций летательных аппаратов. Фотографии образцов реальных 3D КМ и пример построенной модели представлены на рис. 5.

Полученные отдельные сечения поверхности прочности первичного разрушения 3D КМ в выбранных осях представлены на рис. 6(a), 6(b), 6(c).

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

(с) в осях аху ,ayz

Рис. 4. Сечения поверхности прочности 1Р-дисперсно-армированного КМ в осях компонент тензора напряжений

Было показано, что в алгоритме явно могут быть выделены две части: вычисления на хосте и на сопроцессоре.

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

4. Выводы

(1) Использование архитектуры Intel MIC обеспечило существенных прирост производительности при решении поставленной задачи, что позволяет сделать вывод об обоснованности применения

(а) фотографии реальных образцов

(Ь) Модель ячейки пе- (с) поле повреждаемости

риодичности

Рис. 5. Моделирование растяжения по оси ОХ углеродо-углеродных 3Б-армированных КМ в РВС ООБ

таких вычислительных систем для решения прикладных задач рассматриваемого типа.

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

(3) В работе были получены сечения поверхностей прочности ряда моделей КМ. При этом каждое сечение было построено на основе автоматически выбираемых 360-ти простых путей нагру-жения. Сделан вывод о технической возможности реконструкции трехмерных проекций 6-тимерной поверхности прочности с последующей ее визуализацией.

(а) в осях ахх,ауу

(b) в осях ахх, az

(с) в осях аху ,ayz

Рис. 6. Сечения поверхности прочности 3Р-армированно-го КМ в осях компонент тензора напряжений

Благодарности

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

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

[1] N. S. Bakhvalov, G.P. Panasenko, Homogenization: Averaging processes in periodic 'media. Mathematical Problems in the Mechanics of Composite Materials, Mathematics and its Applications, vol. 36, Springer, 1989, URL: http://www.springer.com/mathematics/dynamical+systems/book/978-0-7923-0049-6t3

[2] G. Pavliotis, A. Stuart, Multiscale methods: Averaging and homogenization, Texts in Applied Mathematics, vol. 53, 2008. t 3

[3] Ю. И. Димитриенко, А. И. Кашкаров, А. А. Макашов. «Разработка конечно-элементного метода решения локальных задач теории упругости на ячейке периодичности для композитов с периодической пространственной структурой», Математика в современном мире, ред. Ю. А. Дробышев, 2004, с. 177-191. t 3

[4] Э. Митчел, Р. Уейт. Метод конечных элементов для уравнений с частными производными, Мир, М., 1981. t 3

[5] О. Зенкевич. Метод конечных элементов в технике, Мир, М., 1975. t 3

[6] А. П. Соколов, В. Н. Щетинин, В. М. Макаренков. «Опыт применения теории графов для создания гибких сопровождаемых масштабируемых программных реализаций сложных вычислительных методов», Материалы XIX Международной конференции по вычислительной механике и современным прикладным программным системам, ВМСППС'2015 (24-31 мая 2015 г., Алушта, Крым, Россия), Издательство МАИ, М., 2015, с. 172-174, URL: http://gcad.bmstu.ru/ sa2pdf/gcdcmp_cfa_VMSPPS2015_B_GrphBasedDev.pdf t 4,11,19

[7] A. K. Malmeister, V. P. Tamuj, G. A. Teters. The resistance of polymer and composite materials, Zinatne, Riga, 1980. t 5,7,8

[8] Б. Е. Победря. «Критерии прочности анизотропного материала», Прикладная математика и механика, 1 (1988), с. 141-144. t 5

[9] Дж. Сендецки. Механика композиционных материалов, Мир, М., 1978. t 5

[10] A. V. Ilinykh, V. E. Vildeman. "Modeling of structure and failure processes of granular composites", Computational continuum mechanics, 5:4 (2012), pp. 443-451. t 5

[11] Е. Ю. Макарова, Ю. В. Соколкин. «Нелинейные многоуровневые модели механики деформирования и разрушения композитов», Механика композиционных материалов и конструкций, Тезисы докладов IV-го Всероссийского симпозиума (4-6 декабря 2012, Москва), с. 51. t 5

[12] 1.1. Goldenblat, V. A. Kopnov. "Strength of glass-reinforced plastics in the complex stress state", Polymer Mechanics, 1:2 (1966), pp. 54-59. t 5

[13] Б. Е. Победря. Механика композиционных материалов, МГУ, М., 1984. t 9,10

[14] Ю. И. Димитриенко, А. П. Соколов. «Многомасштабное моделирование упругих композиционных материалов», Математическое моделирование, 24:5 (2012), с. 3-20. t 9

[15] Н. С. Бахвалов. «Осредненные характеристики тел с периодической структурой», Докл. АН СССР, 218:5 (1974), с. 1040. t 10

[16] А. П. Соколов. Математическое моделирование эффективных упругих композитов с многоуровневой иерархической структурой, Ph. D. Thesis, МГТУ им. Н. Э. Баумана, М., 2008. t 10

Рекомендовал к публикации Программный комитет

Четвёртого национального суперкомпьютерного форума НСКФ-2015

Об авторах:

Александр Павлович Соколов

Доцент по кафедрам «Вычислительная математика и математическая физика» и «Системы автоматизированного проектирования», читаемые курсы: «Основы программирования на С++», «Базы данных», «Основы метода конечных элементов», «Системное и прикладное ПО». Веду научно-исследовательскую деятельность в МГТУ им. Н.Э. Баумана по направлению: «Разработка распределенных вычислительных систем анализа свойств композиционных материалов». Автор Распределенной вычислительной системы GCD.

e-mail: [email protected]

Виталий Николаевич Щетинин

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

e-mail: [email protected]

Арсений Сергеевич Сапелкин

Студент кафедры «Вычислительная математика и математическая физика». Область научных интересов: микромеханика композиционных материалов, разработка программных инструментов автоматического построения трехмерных моделей, высокопроизводительные вычисления, разработка распределенных вычислительных систем.

e-mail: [email protected]

Пример ссылки на эту публикацию:

А. П. Соколов, В. Н. Щетинин, А. С. Сапелкин. «Параллельный алгоритм реконструкции поверхности прочности композиционных материалов для архитектуры Intel MIC (Intel Many Integrated Core Architecture)», Программные системы: теория и приложения, 2016, 7:2(29), с. 3-25. URL: http://psta.psiras.ru/read/psta2016_2_3-25.pdf

Alexandr Sokolov, Vitaly Schetinin, Arseniy Sapelkin. Strength surface reconstruction using special 'parallel algorithm based on Intel MIC (Intel Many Integrated Core) architecture.

Abstract. The main objective of this research was to create a parallel software implementation of a numerical method for the reconstruction of strength surface of initial failure of composite material under research. Quadratic tensor polynomial criterion named after Malmeyster-Woo was used. Asymptotic averaging method or homogenization (Nikolay S. Bahvalov, Boris E. Pobedria) and finite element method were used. Software implementation was developed using graph-based technology provided in Distributed computational system GCD using Intel MIC architecture. Computational experiments have been carried out for a series of models of composite materials defined by their reinforcement schemes ("unit cells") to obtain the slice of its strength surface. Three types of composites were considered: 3D-reinforced and ID-reinforced. The results of calculations are presented. (In Russian).

Key words and phrases: micromechanics of composite materials, elastic and strength properties, quadratic failure criterion, homogenization, numerical methods for solving inverse problems, high performance computing, distributed systems, finite element method, the method of asymptotic averaging.

References

[1] N. S. Bakhvalov, G. P. Panasenko, Homogenization: Averaging processes in periodic media. Mathematical Problems in the Mechanics of Composite Materials, Mathematics and its Applications, vol. 36, Springer, 1989, URL: http://www.springer.com/mathematics/dynamical+systems/book/978-0-7923-0049-6

[2] G. Pavliotis, A. Stuart, Multiscale methods: Averaging and homogenization, Texts in Applied Mathematics, vol. 53, 2008.

[3] Yu. I. Dimitriyenko, A. I. Kashkarov, A. A. Makashov. "Razrabotka konechno-elementnogo metoda resheniya lokal'nykh zadach teorii uprugosti na yacheyke periodichnosti dlya kompozitov s periodicheskoy prostranstvennoy strukturoy", Matematika v sovremennom mire, ed. Yu. A. Drobyshev, 2004, pp. 177—191.

[4] E. Mitchel, R. Uyeyt. Metod konechnykh elementov dlya uravneniy s chastnymi proizvodnymi, Mir, M., 1981.

[5] O. Zenkevich. Metod konechnykh elementov v tekhnike, Mir, M., 1975.

[6] A. P. Sokolov, V. N. Shchetinin, V. M. Makarenkov. "Opyt primeneniya teorii grafov dlya sozdaniya gibkikh soprovozhdayemykh masshtabiruyemykh programmnykh realizatsiy slozhnykh vychislitel'nykh metodov", Materialy XIX Mezhdunarodnoy konferentsii po vychislitel'noy mekhanike i sovremennym prikladnym programmnym sistemam, VMSPPS'2015 (24—31 maya 2015 g., Alushta, Krym, Rossiya), Izdatel'stvo MAI, M., 2015, pp. 172-174, URL: http://gcad.bmstu.ru/sa2pdf/gcdcmp_cfa_VMSPPS2015_B_GrphBasedDev.pdf

[7] A. K. Malmeister, V. P. Tamuj, G.A. Teters. The resistance of polymer and composite materials, Zinatne, Riga, 1980.

[8] B.Ye. Pobedrya. "Kriterii prochnosti anizotropnogo materiala", Prikladnaya m,atem,atika i mekhanika, 1 (1988), pp. 141-144.

© A. Sokolov, V. Schetinin, A. Sapelkin, 2016 © Institute of Cytology and Genetics SB RAS, 2016 © Program systems: Theory and Applications, 2016

[9] Dzh. Sendetski. Mekhanika kompozitsionnykh materialov, Mir, M., 1978.

[10] A. V. Ilinykh, V. E. Vildeman. "Modeling of structure and failure processes of granular composites", Computational continuum mechanics, 5:4 (2012), pp. 443-451.

[11] Ye. Yu.Makarova,Yu. V.Sokolkin. "Nelineynyyemnogourovnevyyemodelimekhaniki deformirovaniya i razrusheniya kompozitov", Mekhanika kompozitsionnykh materialov i konstruktsiy, Tezisy dokladov IV-go Vserossiyskogo simpoziuma (4-6 dekabrya 2012, Moskva), pp. 51.

[12] I.I. Goldenblat, V. A. Kopnov. "Strength of glass-reinforced plastics in the complex stress state", Polymer Mechanics, 1:2 (1966), pp. 54-59.

[13] B. Ye. Pobedrya. Mekhanika kompozitsionnykh materialov, MGU, M., 1984.

[14] Yu. I. Dimitriyenko, A. P. Sokolov. "Mnogomasshtabnoye modelirovaniye uprugikh kompozitsionnykh materialov", Matematicheskoye modelirovaniye, 24:5 (2012), pp. 3-20.

[15] N. S. Bakhvalov. "Osrednennyye kharakteristiki tel s periodicheskoy strukturoy", Dokl. ANSSSR, 218:5 (1974), pp. 1040.

[16] A. P. Sokolov. Matematicheskoye modelirovaniye effektivnykh uprugikh kompozitov s mnogourovnevoy iyerarkhicheskoy strukturoy, Ph. D. Thesis, MGTU im. N. E. Baumana, M., 2008.

Sample citation of this publication:

Alexandr Sokolov, Vitaly Schetinin, Arseniy Sapelkin. "Strength surface reconstruction using special parallel algorithm based on Intel MIC (Intel Many Integrated Core) architecture", Program systems: theory and applications, 2016, 7:2(29), pp. 3-25. (In Russian).

URL: http://psta.psiras.ru/read/psta2016_2_3-25.pdf