CC BY
61
УДК 004. 051
ПРИМЕНЕНИЕ ТЕХНОЛОГИИ ПАРАЛЛЕЛЬНЫХ ВЫЧИСЛЕНИЙ ПРИ РАСЧЕТЕ ПЛОСКИХ СТАТИЧЕСКИХ ЗАДАЧ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ
© Нгуен Зуи Тхаи1
Иркутский государственный технический университет,
664074, Россия, г. Иркутск, ул. Лермонтова, 83.
Рассматривается задача расчёта плоских статических задач, в частности расчёта плоских ферм методом конечных элементов. При этом перемещения узлов фермы получаются путем решения системы уравнений равнове-
1 Нгуен Зуи Тхаи, аспирант, тел.: +79246256668, e-mail: duythaipistu@gmail.com Nguyen Duy Thai, Postgraduate, tel.: +79246256668, e-mail: duythaipistu@gmail.com
ВЕСТНИК ИрГТУ №9 (80) 2013
17
сия. Подробно описываются процессы заполнения матрицы жесткости отдельного конечного элемента, матрицы жесткости фермы, матрицы-столбца неизвестных перемещений и матрицы-столбца внешних сил. Предлагается способ применения технологии параллельных вычислений системы MATLAB для повышения эффективности расчёта путем увеличения скорости вычислений.
Ил. 3. Табл. 2. Библиогр. 6 назв.
Ключевые слова: параллельные вычисления; расчет плоских статических задач; расчет плоских ферм; метод конечных элементов; система MATLAB.
PARALLEL COMPUTING TECHNOLOGY APPLICATION IN PLANE STATIC PROBLEMS CALCULATION BY FINITE ELEMENT METHOD Nguyen Duy Thai
Irkutsk State Technical University 83 Lermontov St., Irkutsk, 664074, Russia.
The article examines the problem of plane static problem calculation, in particular the calculation of plane truss by the finite element method. The movements of truss nodes are received by solving the system of equilibrium equations. A detailed description is given to the processes of filling the stiffness matrix of an individual finite element, a truss stiffness matrix, a column matrix of unknown movements and a column matrix of external forces. The author proposes the method of using the parallel computing technology of MATLAB system to improve calculation efficiency by increasing calculation speed.
3 figures. 2 tables. 6 sources.
Key words: parallel computing; calculation of plane static problems; calculation of plane farms; finite element method; MATLAB system.
Введение. В настоящее время метод конечных элементов (МКЭ) представляет собой эффективный численный метод решения инженерных и физических задач. Область его применения простирается от анализа напряжений в конструкциях самолётов или автомобилей до расчёта таких сложных систем, как атомная электростанция. С его помощью рассматривается движение жидкости по трубам, через плотины, в пористых средах, исследуется течение сжимаемого газа, решаются задачи электростатики и смазки, анализируются колебания систем [2]. Кроме этого, с его помощью проводятся расчеты по определению напряженно-деформированного состояния и несущей способности реальных конструкций в самых различных отраслях техники, строительства и машиностроения [1].
В частных случаях машиностроения МКЭ позволяет решать ряд задач, связанных с расчётами стержневых систем, работающих на растяжение-сжатие, прочностью валов, балок, пластин, плоских рам, стержней, испытывающих сложный изгиб, задачи теории упругости, реализации напряженно-деформированного состояния конструкций, переноса тепла и др. [2].
При расчёте стержневых систем, работающих на растяжение, в частности, при расчёте плоских ферм, основным в МКЭ является формирование матрицы жесткости элементов, которая устанавливает связь между узловыми силами элемента и перемещениями узлов. Конечно-элементная сетка такой задачи может иметь большое количество узлов, достигающих иногда одного миллиона. Для решения такой задачи затрачивается большое количество времени. Поэтому возникает проблема ускорения времени расчёта плоских ферм, и этот вопрос можно решать с применением технологии параллельных вычислений.
В данной работе описывается реализация анализа МКЭ для расчёта плоских ферм и применение технологии параллельных вычислений в системе MATLAB для ускорения решения этой задачи.
Постановка задач. Дана плоская ферма, схема которой изображена на рис. 1. Предполагаются следующие параметры фермы: высота фермы на опоре I п. = 2200 мм; высота фермы по середине I т = 3200 мм;
пролет L = 24000 мм; количество панелей по верхнему поясу n = 4; площадь сечения всех стержней
Ak = 1600 мм2 и модуль упругости материала стержней E = 2.06-105 Н/мм2.
Из рис. 1 видно, что ферма состоит из 8 узлов, образующих 13 стержней. Весом стержней фермы пренебрегаем. На узлы 4-8 действуют внешние силы по оси OY и равные 105 Н по модулю. В узле 1 устанавливаются внешние связи по осям OX и OY (заделка), а в узле 3 устанавливается только внешняя связь по оси OY (опора). Задача состоит в определении перемещений узлов и реакций во всех стержнях фермы.
Метод конечных элементов в расчёте плоских ферм. Основным положением в МКЭ являются: формирование матрицы жесткости фермы, которая устанавливает связь между узловыми силами элемента и перемещениями его узлов; решение системы уравнений равновесия, которая имеет следующий вид:
K ■ Z - P = 0, (1)
где K - глобальная матрица жесткости фермы; Z - матрица-столбец неизвестных перемещений; P - матрица-столбец внешних сил.
Рис. 1. Графическая схема фермы
Нахождение матрицы жесткости каждого конечного элемента. При количестве узлов фермы т размеры матрицы жесткости фермы К будут 2т х 2т. Матрица К образуется из всех матриц жесткости отдельных конечных элементов и их количество в данном случае равно 13 матрицам. Обозначим матрицы жесткости конечного элемента К, * = 1,13 в соответствии с номерами стержней фермы.
Для каждого стержня фермы в общей системе координат его матрица жесткости имеет размеры 4x4 и следующий вид:
К =
ЕА
ЕА
■ 0082 0
ЕА
ЕА
■81П0 ■ 0080
■ 81П0^ 0080
ЕА
■ 81П2 0
ЕА
-■ 81П0 ■ 0080
ЕА
ЕА
ЕА
■81П0 ■ 0080
ЕА
I
■81П0 ■ 0080
■ 81П0 ■ 0080
к
ЕА
к
■ 0082 0
ЕА
- ■ 81П 0 ■ 0080
ЕА
■ 81П2 0
ЕА
■81П0 ■ 0080
ЕА
(2)
где Е - модуль упругости материала стержней (Е = 2.06-105 Н/мм2); А - площадь сечения стержня (А = 1600 мм2); 1к - длина стержня; 0 - угол наклона стержня к горизонтальной оси ОХ.
При заданных величинах стержень 1, соединяющий узлы 1 и 2, имеет длину /х = 12000 мм и угол наклона 0 = 0. Запишем матрицу жесткости К стержня (1):
К = 105
2.7467 0 - 2.7467 0
0 0 0 0
- 2.7467 0 2.7467 0
0000
Матрицу К можно записать в виде блоков:
К = 105
2.7467 0 -2.7467 0"
0 0 0 0 " 1,1 1,2"
-2.7467 0 2.7467 0 2,1 2,2
0 0 0 0
Нахождения глобальной матрицы жесткости фермы. Глобальная матрица К, размеры которой равны 16х16, вначале состоит только из нулей. Матрицу К также можно записать в виде блоков:
К =
1,1 2,1 3,1 4,1 5,1 6,1 7,1 8,1
1,2 2,2 3,2 4,2 5,2 6,2 7,2 8,2
1,3 2,3 3,3 4,3 5,3 6,3 7,3 8,3
1,4 2,4 3,4 4,4 5,4 6,4 7,4 8,4
1,5 2,5 3,5 4,5 5,5 6,5 7,5 8,5
1,6 2,6 3,6 4,6 5,6 6,6 7,6 8,6
1,7 2,7 3,7 4,7 5,7 6,7 7,7 8,7
1,8 2,8 3,8 4,8 5,8 6,8 7,8 8,8
= 0
Заполним матрицу К матрицей К путем сложения по элементам соответствующих блоков: (1,1), (1,2), (2,1) и (2,2) матриц К и К и получим новую матрицу К:
К = К + К = 105
467 0 - 2.7467 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
7467 0 2.7467 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Аналогично заполним матрицу К остальными матрицами К., і = 2,13 и в результате получим глобальную матрицу жесткости фермы К (количество знаков после запятой принимает 2):
К = 10'
0.69 0.19 -0.27 0 0 0 0 0 -0.42 -0.19 0 0 0 0 0 0
0.19 1.58 0 0 0 0 0 -1.50 -0.19 -0.08 0 0 0 0 0 0
- 0.27 0 1.38 0 -0.27 0 0 0 -0.41 0.19 0 0 -0.42 -0.1875 0 0
0 0 0 1.20 0 0 0 0 0.19 -0.08 0 -1.03 -0.19 -0.0844 0 0
0 0 -0.27 0 0.69 -0.19 0 0 0 0 0 0 -0.42 0.1875 0 0
0 0 0 0 -0.19 1.58 0 0 0 0 0 0 0.19 -0.0844 0 -1.50
0 0 0 0 0 0 0.54 -0.54 -0.54 -0.04 0 0 0 0 0 0
0 -1.50 0 0 0 0 0.04 -0.04 -0.04 0 0 0 0 0 0 0
- 0.42 -0.19 - 0.42 0.19 0 0 -0.54 1.92 1.92 0.09 -0.54 -0.04 0 0 0 0
-0.19 -0.08 0.19 -0.09 0 0 -0.04 0.09 0.09 0.18 -0.04 0 0 0 0 0
0 0 0 0 0 0 0 -0.54 -0.54 -0.04 1.09 0 -0.54 0.0453 0 0
0 0 0 -1.04 0 0 0 -0.04 -0.04 0 0 1.04 0.04 -0.0038 0 0
0 0 - 0.42 -0.19 - 0.42 0.19 0 0 0 0 -0.54 0.04 1.92 -0.0906 -0.54 0.04
0 0 -0.19 -0.08 0.19 -0.08 0 0 0 0 0.04 0 -0.09 0.1763 0.04 0
0 0 0 0 0 0 0 0 0 0 0 0 -0.54 0.0453 0.54 -0.04
0 0 0 0 0 -1.50 0 0 0 0 0 0 0.04 -0.0038 -0.04 1.50
Матрица-столбец неизвестных перемещений и внешних сил. Матрица-столбец неизвестных перемещений является последовательным перечислением перемещений вдоль осей ОХ и OY каждого узла фермы:
гг=[и г, и, V из V и4 V. и, V, и V и7 у, и, у,],
где и и V - перемещения нго узла вдоль осей ОХ и OY.
Внешние нагрузки приложены к узлам 4-8 против направления оси OY. По всем другим направлениям и в других узлах внешние силы равны нулю. Матрица-столбец внешних сил выглядит как последовательное перечисление сил вдоль осей ОХ и OY каждого узла фермы:
Рт =[Х, Г Х2 7 X, Г X4 Г4 X, 7, X Г X, У7 X, У, ],
где X- и у - силы нго узла вдоль осей ОХ и OY.
Рт = [0 0 0 0 0 0 0 -104 0 -104 0 -104 0 -104 0 —104].
Нахождение перемещений узлов с учётом внешних связей. Перепишем систему уравнений для определения перемещений узлов фермы в виде (1). Заметим, что перемещение вдоль оси 0Х узла 1 равно нулю, так как по этому направлению в этом узле установлена линейная связь. Это означает, что соответствующие строка и столбец могут быть вычеркнуты. По такому же принципу, учитывая отсутствие внешних связей в узлах 2 и 6, вычеркнем из системы уравнений (1) строки и столбцы с номерами 2 и 6. Получим окончательную систему уравнений:
105 ■
1.38 0 -0.27 0 0 -0.42 0.19 0 0
0 0.20 0 0 0 0.19 -0.08 0 -1.03
-0.27 0 0.69 0 0 0 0 0 0
0 0 0 0.54 - 0.54 -0.54 -0.04 0 0
0 0 0 0.04 -0.04 -0.04 0 0 0
-0.42 0.19 0 - 0.54 1.92 1.92 0.09 -0.54 -0.04
0.19 -0.08 0 - 0.04 0.09 0.09 0.18 -0.04 0
0 0 0 0 - 0.54 -0.54 -0.04 1.09 0
0 -1.03 0 0 -0.04 -0.04 0 0 1.04
-0.42 -0.19 -0.42 0 0 0 0 -0.54 0.04
-0.19 -0.08 0.19 0 0 0 0 0.04 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
- 0.42 - 0.19 0
0
-0.42 0.19
-0.,4 0.04 0.04 0
0
0
1.92 - 0.09 - 0.,4 - 0.09 0.18 0.04
-0.,4 0.04 0.04 0
0
0.04
0
0.,4 -0.04 -0.04 1.,0
и2 ' 0 '0'
У2 0 0
и3 0 0
и4 0 0
У4 -104 0
и5 0 0
- -104 = 0
иь 0 0
-104 0
и7 0 0
У7 -104 0
и8 0 0
У8 -104 0
Решение системы уравнений, например, методом исключения Гаусса дает следующие искомые параметры матрицы-столбца неизвестных перемещений Z : и2 = 1,2136, У2 =-7,2030, из = 2,4272, и4 = 1,3056,
V =-0,0667, и = 1,7774, V =-5,7279, иб = 1,2136, У6 =-7,2394, и7 = 0,6498, У7 =-5,7279,
и8 = 1,1215 , V =-0,0667 .
Сведем величины в табл.1.
Перемещения узлов фермы
Таблица 1
Номер узла Перемещение Номер узла Перемещение
ОХ, мм OY, мм ОХ, мм OY, мм
1 0 0 5 1,7774 -5,7279
2 1,2136 -7,2030 6 1,2136 -7,2394
3 2,4272 0 7 0,6498 -5,7279
4 13056 -0,0667 8 1,1215 -0,0667
0
0
0
0
0
0
0
0
0
0
0
0
0
Вычисление усилий в стержнях. Для вычисления усилий в стержне нужно использовать следующую последовательность действий.
1. Определение новой длины стержня:
/, — л/ (X — X ) + (у — у ) ,
к —новая X \ кн нч / кн у нч / '
где хт , _укн - координаты начального узла (любого из двух) стержня после деформирования; хш , - коор-
динаты конечного узла (отличного от начального) стержня после деформирования.
2. Определение нового угла наклона стержня: оновый — arctan[( ун — унч)/(хкн — хт )].
3. Определение разности углов наклона по модулю: |До| — |о — ожвый|.
4. Нахождение проекции новой длины стержня: 4_ироеи<мм — /
Д/„
к—новая
• cos(До).
5. Нахождение деформации стержни: є — ■
/
к — 4
к—проекции к
кк
6. Вычисление усилия в стержне: N — є • Е • ^.
Результат вычисления усилий в стержнях показан в табл. 2 и на рис. 2, где изображена графическая схема фермы после деформации.
Таблица 2
Усилия в стержнях
Номер стержня 2 Усилие, Н/мм Номер стержня 2 Усилие, Н/мм
1 33333,3333 8 4569,1067
2 33333,3333 9 4569,1067
3 0 10 -36552,8537
4 -37629,9831 11 -10000
5 -37629,9831 12 -3570
6 0 13 -10000
7 -36552,8537
Рис. 2. Графическая схема фермы после деформации
При расчёте плоских ферм в инженерных расчётах с помощью программных систем, например, программной системы МА^АВ, предлагаются последовательный и параллельный варианты алгоритма.
Последовательный вариант алгоритма расчёта плоских ферм. Последовательный вариант алгоритма расчёта плоских ферм состоит из следующих шагов:
- шаг 1: ввод начальных значений переменных Ноп, Н , Ь , п , Е, А;
- шаг 2: назначение внешних сил и ввод их величин;
- шаг 3: назначение внешних связей;
- шаг 4: определение количества узлов фермы и их координат;
- шаг 5: определение количества стержней фермы и номеров конечных узлов;
- шаг 6: вычисление длин и углов наклона стержней;
- шаг 7: заполнение глобальной матрицы жесткости фермы матрицами жесткости стержней;
- шаг 8: вычеркивание строк и столбцов глобальной матрицы жесткости фермы, матрицы неизвестных пере-
мещений и матрицы внешних сил с учётом внешних связей;
- шаг 9: вычисление перемещений узлов путем решения системы (1);
- шаг 10: вычисление новых длин и углов наклона стержней;
- шаг 11: вычисление усилий во всех стержнях;
- шаг 12: вывод результатов.
Параллельный вариант алгоритма расчёта плоских ферм. Результаты тестирования последовательного варианта алгоритма расчёта плоских ферм показывают, что большая часть времени тратится на заполнение матриц жесткости отдельных стержней (шаг 7) и решение системы линейных алгебраических уравнений (СЛАУ) для вычисления перемещений узлов (шаг 9). Поэтому предлагается выполнять эти шаги одновременно на нескольких исполнителях (ядрах процессора), то есть с помощью технологии параллельных вычислений. Была выбрана система MATLAB, так как она поддерживает технологию параллельных вычислений с помощью двух пакетов расширения Parallel Computing Toolbox и MATLAB Distributed Computing Server. Технология параллельных вычислений в системе MATLAB была описана [3-6].
Параллельное заполнение матрицы жесткости фермы матрицами жесткости стержней. Для заполнения матрицы жесткости фермы нужно найти матрицы жесткости стержней и добавить их к матрице жесткости фермы. Подход, лежащий в основе параллельного заполнения матрицы жесткости фермы, заключается в том, что на начальной стадии (инициализации) все стержни (их количество, номера конечных узлов) расположены в рабочем пространстве первого исполнителя (исполнитель с индексом 1). Сначала первый исполнитель рассылает всем исполнителям, включая и себя, равные количества стержней. На каждом исполнителе вычисляются матрицы жесткости своих стержней, которые размещаются в своей глобальной матрице жесткости. После этого все остальные исполнители, кроме первого, отправляют свою глобальную матрицу жесткости первому исполнителю. В результате в рабочем пространстве первого исполнителя будет храниться глобальная матрица жесткости фермы.
Параллельные методы решения системы СЛАУ. Как известно из курса линейной алгебры, популярными методами решения СЛАУ AX = B являются методы исключения Гаусса, простой итерации, Крамера и др. Подход, лежащий в основе параллельного метода решения СЛАУ методом исключения Гаусса, заключается в том, что на начальной стадии (инициализации) вся матрица A расположена в рабочем пространстве первого исполнителя. На первом шаге первый исполнитель отправляет слои матрицы A всем остальным исполнителям. После отправки первый исполнитель нормирует первый элемент первой строки своей подматрицы и сразу же отправляет её всем остальным исполнителям. Остальные исполнители, получив данную строку, начинают исключать неизвестные первые элементы из строк своей подматрицы. К тому моменту, когда первый исполнитель завершил приведение своей подматрицы к верхнетреугольному виду, остальные исполнители также завершают получение от него соответствующих строк. После этого все исполнители отправляют оставшиеся ненулевые подматрицы первому исполнителю. В итоге после первого шага в рабочем пространстве первого исполнителя хранится верхнетреугольная подматрица и квадратная ненулевая подматрица (полученная из подматриц, присланных остальными исполнителями). Следующий шаг заключается в том, что первый исполнитель рассылает соответствующие подматрицы ненулевой квадратной матрицы всем исполнителям, включая и себя. Алгоритм продолжает работать до тех пор, пока размер отправляемого слоя квадратной матрицы не станет равным единице.
Рис. 3. График зависимости временных затрат решения задачи от количества панелей фермы
В результате, после последнего шага алгоритма в рабочем пространстве первого исполнителя будет храниться верхнетреугольная матрица U с единичной диагональю. Именно эта матрица и будет считаться результатом работы параллельного алгоритма исключения Гаусса.
Подобные описания подходов, лежащих в основе параллельных методов решения СЛАУ, описаны в [3, 4]. Экспериментальная оценка эффективности. Критерий эффективности представляет собой разность времен At, затрачиваемых для расчета в последовательном и параллельном вариантах алгоритма.
Для проведения тестирования алгоритмов был использован локальный планировщик, запущенный на ноутбуке Alienware со следующей конфигурацией:
- процессор Intel Core i7 - 2720 - 2.2 Гц с 8 ядрами в каждом;
- оперативная память DDR3 с размером 8Гб и частотой 1600 МГц;
- объём жесткого диска 750 Гб;
- операционная система Window 7 64-bit.
На рис. 3 изображены графики зависимости временных затрат решения задачи при использовании одного и двух исполнителей и изменении количества панелей фермы от 50 до 500.
Заключение. По результатам тестирования можно сделать вывод, что при количестве панелей больше 250-ти, время, затрачиваемое на расчёт в параллельном режиме, меньше, чем в последовательном при количестве панелей n = 500 и количестве исполнителей p = 2, At = 14.53 с. При увеличении количества исполнителей от 2 до 8 время, затрачиваемое на расчёт, увеличивается. Это можно объяснить тем, что время, затрачиваемое на обмен данными между исполнителями, превышает выигрыш от быстродействия вычислений. На практике, при задаче, конечно-элементная сетка которой может иметь большое количество узлов, достигающих иногда одного миллиона, применение технологии параллельных вычислений позволяет ускорить время расчёта, т.е. эффективность имеет место.
Библиографический список
1. Куканов Н.И., Черный А.Н. Расчет фермы методом конечных элементов: метод. указания. Ульяновск: Изд-во УлГТУ, 2005. 28 с.
2. Овраченко В.А. Расчет задач машиностроения методом конечных элементов: учеб. пособие. Краматорск: ДГМА, 2004. 128 с.
3. Нгуен Зуи Тхаи., Сосинская С.С. Разработка лабораторного практикума по численным методам в параллельном режиме // Фундаментальные и прикладные исследования, разработка и применение высоких технологий в промышленности и экономике: сб. стат. СПб.: Изд-во Политехнического ун-та. 2012. Т. 1. С. 115-120.
4. Оленёв Н.Н. Основы параллельного программирования в системе MPI. М.: ВЦ РАН, 2005. 80 с.
5. Parallel Computing Toolbox User's Guide [Электронный ресурс]. Режим доступа: http://www.mathworks.com/help/pdf_doc/distcomp/distcomp.pdf
6. MATLAB Distributed Computing Server System Administrator's Guide [Электронный ресурс]. Режим доступа: http://www.mathworks.com/help/pdf_doc/mdce/mdce.pdf
