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

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

Рисунок 1. Заданная расчетная схема рамы.

1. Разработка схемы дискретизации.

Разработка схемы дискретизации включает в себя следующие действия:

1.1. Обозначение ГСК (глобальной системы координат) XYZ. Начало координат располагаем таким образом, чтобы в дальнейшем все координаты узлов рамы были положительными (рис.2). Ось X направлена вправо, Y - вверх, Z - смотрит на нас.

Рисунок 2. Обозначение глобальной системы координат.

1.2. Нумерация узлов рамы. Узлы рамы обычно нумеруют слева направо и снизу вверх. В рассматриваемом примере, двигаясь слева направо, нумеруем узлы 1 и 2. Затем двигаемся снизу вверх (когда узлы рамы расположены на одной линии по вертикали), нумеруем сначала нижний узел 3, затем – верхние узлы 4 и 5. Двигаясь далее слева направо, нумеруем последний узел 6 (рис.3). В рамах с большим количеством узлов нужно стремиться к тому, чтобы разность между номерами узлов одного конечного элемента была минимальной.

Рисунок 3. Рама с пронумерованными узлами.


 

1.3. Нумерация конечных элементов соответствует нумерации узлов рамы. Первым нумеруем конечный элемент, который берет свое начало в точке 1. При этом конечный номер этого элемента должен быть минимальным. То есть, первый элемент 1-2.

Затем нумеруем конечные элементы, имеющие начало в точке 2. При этом, сперва нумеруем элемент, конечная точка которого 3, и только после этого – элемент, конечная точка которого 4. Таким образом, получаем второй элемент 2-3 и третий элемент 2-4.

 Далее нумеруем конечный элемент, берущий начало в точке 4, и имеющий конец в точке 5. Это четверной элемент 4-5.

Соответственно, последним нумеруем пятый элемент 5-6 (рис.4).

Рисунок 4. Рама с пронумерованными КЭ.

1.4. Обозначение локальной системы координат для каждого КЭ - xyz. Начало каждой из локальных систем координат находится в начальном узле рассматриваемого конечного элемента. Ось x направлена вдоль КЭ, ось y - против часовой стрелки под углом 90° к оси x, ось z (на рис. 5 условно не показана) - смотрит на нас.

Рисунок 5. ЛСК, обозначенные для каждого КЭ.

1.5. Обозначение направляющих углов φ для каждого КЭ. Этот угол определяется от оси X  глобальной системы координат против часовой стрелки до совмещения с осью x локальной системы координат каждого конечного элемента (рис.6).

Рисунок 6. Направляющие углы, обозначенные для каждого КЭ.


 

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

Рисунок 7. Схема дискретизации

Таблица 1

Координаты узлов дискретизации

№ узла 1 2 3 4 5 6
х, а 0 2 2 2 2 4
у, а 2 2 0 3 4 4

 

Определение вспомогательных величин проводят по формулам аналитической геометрии:

i - номер начального узла КЭ; j - номер конечного узла КЭ; ;

; ; .

Заметим, что в рассматриваемом примере определение длин конечных элементов, а также косинусов (синусов) направляющих углов может быть найдено и без применения формул.

Так, например, длины элементов легко определяются непосредственно по рис. 7, а значение косинуса (синуса) углов 0⁰ и 90⁰ должно быть хорошо известно из курса тригонометрии:

  

Вычисления заносят в таблицу 2.

Таблица 2

Обработка узлов дискретизации

№ КЭ i j lk, a sin jk cos jk
1 1 2 2 0 1
2 2 3 2 -1 0
3 2 4 1 1 0
4 4 5 1 1 0
5 5 6 2 0 1

 


 

Матрицы преобразования для балочного конечного элемента имеют вид:

Следовательно, для поставленной задачи, пользуясь данными таблицы 2, имеем:

3. Формирование матриц жесткости в глобальной системе координат.

В зависимости от используемого типа балочного конечного элемента, матрица жесткости в локальной системе координат записывается следующим образом (EJ – жесткость стержней на изгиб, в рассматриваемом примере жесткость горизонтальных элементов принята в 2 раза больше, чем у вертикальных):

 


 

Разбиваем ЗРС на балочные конечные элементы двух типов (рис. 8):

Рисунок 8. Разбиение ЗРС на балочные КЭ двух типов

Пользуясь данными таблицы 2, учитывая тип КЭ, получаем:

4. Вычисление реакций на внутрипролетную нагрузку, приложенную к КЭ.

4.1. Из ЗРС видно, что ко второму КЭ приложена распределенная нагрузка. Обратившись к табличным эпюрам метода перемещений, можно определить значеня элементо вектора реакций в виде:

4.2. Ко второму КЭ также приложен момент. Обратившись к табличным эпюрам метода перемещений, можно определить значения элементов вектора реакций в виде:

Сложив значения эпюры от действия распределенной нагрузки и эпюры от действия момента, получаем результирующую эпюру моментов для второго КЭ, которая будет иметь вид:

 


 

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

Далее по построенным эпюрам для каждого КЭ составляем матрицы реакций на внутрипролетную нагрузку в локальной системе координат. При этом необходимо учитывать, что значения реакций записываются в матрицы с учетом математического правила знаков (рис.9).

Рисунок 9. Схема, иллюстрирующая математическое правило знаков.

В том случае, если реакция в соответствующей точке рассматриваемого КЭ оказывается направлена в ту же сторону, что и на рисунке 9, то значение данной реакции берется с положительным знаком. В противном случае значение данной реакции следует записывать в матрицу с отрицательным знаком.


 

5. Матрицы жесткости из локальной системы координат преобразовываются в глобальную систему координат по формуле:

6. Матрицы реакций на внутрипролетную нагрузку в глобальной системе координат определяются по формуле:

7. Формирование матрицы жесткости для ансамбля КЭ в ГСК.

Обозначим перемещения узлов рамы в ГСК (рис.10). Для этого в каждом узле рамы обозначим по три степени свободы – горизонтальную, вертикальную и угловую. Пока значения этих перемещений неизвестны, они принимаются положительными, то есть направленными вдоль положительных полуосей ГСК, поворот – против часовой стрелки.

Рисунок 10. Обозначение степеней свободы в каждом узле рамы

Так как на раму наложены внешние связи, прикрепляющие ее к плоскости, то некоторые из этих перемещений отсутствуют. В точке 1 запрещено горизонтальное перемещение, в точке 3 - вертикальное, в точке 6 запрещены три перемещения – горизонтальное, вертикальное и угловое.

Из ЗРС видно, что горизонтальных перемещений не будет также в точках 2, 5, вертикальных – в точках 2, 4, 5.

Для упрощения дальнейших расчетов, нужно стремиться к тому, чтобы количество перемещений было минимальным. Анализируя рис.8, а также учитывая типы используемых балочных конечных элементов, понимаем, что угловые перемещения будут отсутствовать в точках 1, 3, 4.

Кроме этого, перемещениями, возникающими вследствие растяжения-сжатия КЭ, пренебрегаем. Таким образом, получаем следующие перемещения рамы в ГСК (рис. 11): первое - вертикальное в точке 1, второе – угловое в точке 2, третье – горизонтальное в точке 3, четвертое - горизонтальное в точке 4 и пятое – угловое в точке 5. Следует заметить, что в одной точке сначала нумеруется горизонтальное, потом вертикальное, а затем угловое перемещение.


 

Рисунок 11. Векторы перемещений в глобальной системе координат.

Таблица 3

Матрица индексов ансамбля КЭ

 

 

№ КЭ

d1 d2 d3 d4 d5 d6
1 0 1 0 0 0 2
2 0 0 2 3 0 0
3 0 0 2 4 0 0
4 4 0 0 0 0 5
5 0 0 5 0 0 0

 

Таблица индексов заполняется следующим образом:

 








Первая строка

- d1 – это горизонтальное перемещение начальной точки первого элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d2 – это вертикальное перемещение начальной точки первого элемента, так как это перемещение имеет номер 1, то на этой позиции ставим 1;

- d3 – это угловое перемещение начальной точки первого элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d4 – это горизонтальное перемещение конечной точки первого элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d5 – это вертикальное перемещение конечной точки первого элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d6 – это угловое перемещение конечной точки первого элемента, так как это перемещение имеет номер 2, то на этой позиции ставим 2;

 

Вторая строка

- d1 – это горизонтальное перемещение начальной точки второго элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d2 – это вертикальное перемещение начальной точки второго элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d3 – это угловое перемещение начальной точки второго элемента, так как это перемещение имеет номер 2, то на этой позиции ставим 2;

- d4 – это горизонтальное перемещение конечной точки второго элемента, так как это перемещение имеет номер 3, то на этой позиции ставим 3;

- d5 – это вертикальное перемещение конечной точки второго элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d6 – это угловое перемещение конечной точки второго элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

Третья строка

- d1 – это горизонтальное перемещение начальной точки третьего элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d2 – это вертикальное перемещение начальной точки третьего элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d3 – это угловое перемещение начальной точки третьего элемента, так как это перемещение имеет номер 2, то на этой позиции ставим 2;

- d4 – это горизонтальное перемещение конечной точки третьего элемента, так как это перемещение имеет номер 4, то на этой позиции ставим 4;

- d5 – это вертикальное перемещение конечной точки третьего элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d6 – это угловое перемещение конечной точки третьего элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;



Четвертая строка

- d1 – это горизонтальное перемещение начальной точки четвертого элемента, так как это перемещение имеет номер 4, то на этой позиции ставим 4;

- d2 – это вертикальное перемещение начальной точки четвертого элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d3 – это угловое перемещение начальной точки четвертого элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d4 – это горизонтальное перемещение конечной точки четвертого элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d5 – это вертикальное перемещение конечной точки четвертого элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d6 – это угловое перемещение конечной точки четвертого элемента, так как это перемещение имеет номер 5, то на этой позиции ставим 5;

Пятая строка

- d1 – это горизонтальное перемещение начальной точки пятого элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d2 – это вертикальное перемещение начальной точки пятого элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d3 – это угловое перемещение начальной точки пятого элемента, так как это перемещение имеет номер 5, то на этой позиции ставим 5;

- d4 – это горизонтальное перемещение конечной точки пятого элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d5 – это вертикальное перемещение конечной точки пятого элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

- d6 – это угловое перемещение конечной точки пятого элемента, так как это перемещение отсутствует, то на этой позиции ставим 0;

 

Формирование матрицы жесткости.

Матрица жесткости ансамбля конечных элементов имеет размер 5х5 (по числу независимых перемещений в матрице индексов):

.

Элементы матрицы вычисляются следующим образом:

7.1. В матрице [K] для данной задачи получилось 25 элементов. Но с учетом равенства элементов, стоящих на побочных диагонялаях K12 = K21, K13 = K31, K14 = K41, K15 = K51, K23 = K32, K24 = K42, K25 = K52, K34 = K43, K35 = K53, K45 = K54, независимых элементов получается 15. Это элементы, стоящие на главной диагонали K11, K22, K33, K44, K55 и элементы, стоящие на побочных диагоналях K12, K13, K14, K15, K23, K24, K25, K34, K35, K45.

7.2. Элементы матрицы [K] формируются в соответствии с таблицей 3.

Найдем значение элемента K11, стоящего на главной диагонали. Для этого в таблице 3 (матрице индексов) выделим позиции, где стоит цифра 1. В таблице такое место одно: в первой строке, значит свой вклад в определение этого коэффициента даст первый конечный элемент. При этом, в первой строке единица стоит на позиции d2.

Таблица 3

Матрица индексов ансамбля КЭ

№ КЭ d1 d2 d3 d4 d5 d6
1 0 1 0 0 0 2
2 0 0 2 3 0 0
3 0 0 2 4 0 0
4 4 0 0 0 0 5
5 0 0 5 0 0 0

Поэтому элемент K11 вычисляется следующим образом:

Верхний индекс «Г» говорит о том, что мы должны взять соответствующий элемент из матрицы жесткости для отдельного конечного элемента в глобальной системе координат (они вычислялись выше).

Первый индекс у элемента  говорит о том, что единица стоит в первой строке, то есть соответствует первому конечному элементу. Индексы после запятой соответствуют номеру перемещения d2, так как единица стоит во втором столбце матрицы индексов.

Далее этот элемент берется из соответствующей матрицы жесткости первого конечного элемента (индексы после запятой указывают соответственно номер строки и номер стоблца в соответствующей матрице). Не забудьте учесть множитель, который выносится за знак матрицы!

Поэтому получаем:

Таким образом, оконачательно для элемента K11 получаем:

Далее найдем значение элемента K22, также стоящего на главной диагонали. Для этого в таблице 3 (матрице индексов) выделим все позиции, где стоит цифра 2. В таблице таких мест три: в первой, второй и третьей строках, значит свой вклад в определение этого коэффициента дадут конечные элементы первый, второй и третий. При этом, в первой строке двойка стоит на позиции d6, во второй и третьей строках – на позиции d3.

Таблица 3

Матрица индексов ансамбля КЭ

№ КЭ d1 d2 d3 d4 d5 d6
1 0 1 0 0 0 2
2 0 0 2 3 0 0
3 0 0 2 4 0 0
4 4 0 0 0 0 5
5 0 0 5 0 0 0

 

Поэтому элемент K22 вычисляется следующим образом:

Первый индекс у элемента  говорит о том, что двойка стоит в первой строке, то есть соответствует первому конечному элементу. Индексы после запятой соответствуют номеру перемещения d6, так как двойка стоит в шестом столбце матрицы индексов.

Первый индекс у элемента  говорит о том, что следующая двойка стоит во второй строке, то есть соответствует второму конечному элементу. Индексы после запятой соответствуют номеру перемещения d3, так как двойка стоит в третьем столбце матрицы индексов.

Первый индекс у элемента  говорит о том, что ещё одна двойка стоит в третьей строке, то есть соответствует третьему конечному элементу. Индексы после запятой соответствуют номеру перемещения d3, так как двойка стоит в третьем столбце матрицы индексов.

Далее каждый из этих элементов берется из соответствующей матрицы жесткости первого, второго и третьего конечных элементов (индексы после запятой указывают соответственно номер строки и номер стоблца в соответствующей матрице). Не забудьте учесть множитель, который выносится за знак матрицы!

Поэтому получаем:

Таким образом, оконачательно для элемента K22 получаем:


 

Аналогично вычисляются остальные элементы, стоящие на главной диагонали матрицы:

7.3. Найдем значение элемента K12, не принадлежащего главной диагонали матрицы. Для этого в таблице 3 (матрице индексов) выделим позиции, где стоят цифры 1 и 2 в одной строке!  В таблице такое место одно: в первой строке, значит свой вклад в определение этого коэффициента даст первый конечный элемент. При этом, в первой строке единица стоит на позиции d2, а двойка – на позиции d6.

Таблица 3

Матрица индексов ансамбля КЭ

№ КЭ d1 d2 d3 d4 d5 d6
1 0 1 0 0 0 2
2 0 0 2 3 0 0
3 0 0 2 4 0 0
4 4 0 0 0 0 5
5 0 0 5 0 0 0

 

Поэтому элемент K12 вычисляется следующим образом:

Первый индекс у элемента  говорит о том, что пара цифр 1 и 2 стоит в первой строке, то есть соответствует первому конечному элементу. Индексы после запятой говорят о том, что цифра 1 соответствует перемещению d2, так как единица стоит во втором столбце матрицы индексов, а цифра 2 соответствует перемещению d6, так как двойка стоит в шестом столбце матрицы индексов.

Далее этот элемент берется из соответствующей матрицы жесткости первого конечного элемента (индексы после запятой указывают соответственно номер строки и номер стоблца в соответствующей матрице). Не забудьте учесть множитель, который выносится за знак матрицы!

Поэтому получаем:

Таким образом, оконачательно для элемента K12 получаем:

Аналогично вычисляются остальные коэффициенты, не принадлежащие главной диагонали матрицы:

 

8. Формирование вектора реакций на внутрипролетную нагрузку.

Количество составляющих вектора реакций на внутрипролетную нагрузку определяется количеством перемещений в глобальной системе координат. В рассматриваемом примере количество перемещений равно пяти, следовательно, вектор-столбец будет иметь пять строк.

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

Рисунок 12. Схема для формирования вектора реакций на внутрипролетную нагрузку.

Значение реакций, направления которых совпадают с направлением перемещения в соответствующей точке, возьмем с положительным знаком, если же направление реакции и перемещения в рассматриваемой точке не совпадают, значение этой реакции берем с отрицательным знаком:

Разрешающие уравнения метода.

Для решения задачи методом конечных элементов с использованием балочных конечных элементов разрешающая система уравнений имеет вид:

 

 - вектор внешних нагрузок. Сила 3qa оказывается на первой позиции в матрице {P} и берется с отрицательным знаком, поскольку ее направление противоположно направлению перемещения  (рис. 13).

Рисунок 13. Составление вектора внешних нагрузок.

 

Получаем систему линейных алгебраических уравнений метода конечных элементов:

9. Решение СЛАУ.

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

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

Перепишем систему уравнений в матричном виде и решим его методом Гаусса:

1-ую строку делим на 0.75

от 2 строки отнимаем 1 строку, умноженную на 1.5

2-ую строку делим на 4.5

от 1 строки отнимаем 2 строку, умноженную на 2; к 3 строке добавляем 2 строку, умноженную на 0.75; от 4 строки отнимаем 2 строку, умноженную на 3

3-ую строку делим на 0.25

от 1 строки отнимаем 3 строку, умноженную на 13; к 2 строке добавляем 3 строку, умноженную на 16; от 4 строки отнимаем 3 строку, умноженную на 0.5

4-ую строку делим на 3

к 1 строке добавляем 4 строку, умноженную на 2; от 2 строки отнимаем 4 строку, умноженную на 1; от 3 строки отнимаем 4 строку, умноженную на 2; к 5 строке добавляем 4 строку, умноженную на 3

5-ую строку делим на 4

к 1 строке добавляем 5 строку, умноженную на 2; от 2 строки отнимаем 5 строку, умноженную на 1; от 3 строки отнимаем 5 строку, умноженную на 2; к 4 строке добавляем 5 строку, умноженную на 1

Сделаем проверку. Подставим полученное решение в уравнения из системы и выполним вычисления:

Проверка выполнена успешно.

 

Таким образом, в результате решения получаем следующие перемещения:

10. Формирование матрицы узловых перемещений КЭ в ЛСК.

Так как перемещения в ЛСК связаны с перемещениями в ГСК, то потребуются составляющие векторов узловых перемещений в ГСК для каждого КЭ. Для дальшейших рассуждений необходимо использовать решение системы уравнений МКЭ с учетом данных табилцы индексов (табл.3).

Проанализируем, как будут перемещаться узловые точки первого КЭ.

Рассмотрим узел 1. Как было отмечено выше, в данном узле возможно вертикальное перемещение Δ1 (см. рис.11). Поскольку в результате решения СЛАУ значение этого перемещения получилось отрицательным, узел 1 будет перемещаться в направлении, противоположном изначально принятому на рис. 11 на величину 7,1667а.

Далее рассмотрим узел 2. Как было отмечено выше, в данном узле имеется угловое перемещение Δ2 (см. рис.11). Поскольку в результате решения СЛАУ значение этого перемещения получилось положительным, узел 2 будет поворачиваться против хода часовой стрелки, то есть в направлении, изначально принятом на рис. 11 на величину 1,5833а.

После этого на примере первого КЭ покажем, по какому принципу составляется матрица узловых перемещинй КЭ в ГСК.

Первое значение соответствует горизонтальному перемещению начальной точки первого элемента, так как это перемещение отсутствует, то на этой позиции стоит 0;

Второе значение соответствует вертикальному перемещению Δ1 начальной точки первого элемента, так как это перемещение в результате решения СЛАУ получилось равным -7,1667а, то во второй строке записано -7,1667а;

Третье значение соответствует угловому перемещению начальной точки первого элемента, так как это перемещение отсутствует, то на этой позиции стоит 0;

Четвертое значение соответствует горизонтальному перемещению конечной точки первого элемента, так как это перемещение отсутствует, то на этой позиции стоит 0;

Пятое значение соответствует вертикальному перемещению конечной точки первого элемента, так как это перемещение отсутствует, то на этой позиции стоит 0;

Шестое значение соответствует угловому перемещению Δ2 конечной точки первого элемента, так как это перемещение в результате решения СЛАУ получилось равным 1,5833, то во второй строке записано 1,5833;

 

Рассуждая аналогично, составляем матрицы узловых перемещинй остальных КЭ в ГСК:


 

Матрицы узловых перемещений конечных элементов в ЛСК получают по формуле:

 

 

 

 

 


 

11. Вычисление усилий

Усилия вычисляются по формуле:

 


 

Деформированное состояние ЗРС, отвечающее полученному решению, приведено на рис. 14.

Рисунок 14. Схема возможных перемещений.

Форма записи значения усилий КЭ с номером k в точках i,j имеет вид:

и позволяет по полученным значениям осуществить построение эпюр M и Q, отложив величины поперечной силы и момента в соответствующих точках i, j конечного элемента. После построения эпюр M, Q осуществляется построение эпюры N.

Результат построения эпюр показан на рис. 15.

 

Рисунок 15. Эпюры М, Q, N.


 

12. Контроль проведенных построений осуществляется проверкой условий равновесия узлов ЗРС и ее произвольной части.

 

 

Рисунок 16. Равновесие конечных элементов (сечение проведены бесконечно близко к узлам).






Узел 1

;

;

Узел 2

;

;

;

Узел 3

;

;

Узел 5

;

;

;

Узел 6

;

;

;


 

Дополнительный контроль проведенных построений может быть осуществлен в программно-вычислительном комплексе SCAD:

Рисунок 17. Эпюра M, qa2.

Рисунок 18. Эпюра Q, qa.

 

Рисунок 19. Эпюра N, qa.

Таблица 4

Сравнение значений усилия ручного расчета и расчета, выполненного
в программе SCAD

№ КЭ

SCAD

Ручной расчет

Расхождение, %

М Q N М Q N М Q N
1 6,0 -3,0 -2,0 6,0 -3,0 -2,0 0 0 0
2 5,0 -2,0 -1,5 5,0 -2,0 -1,5 0 0 0
3 1,0 -2,0 1,5 1,0 -2,0 1,5 0 0 0
4 2,0 -2,0 1,5 2,0 -2,0 1,5 0 0 0
5 2,0 -1,5 2,0 2,0 -1,5 2,0 0 0 0

 



Вывод

Проверка полученных результатов показывает, что условия равновесия выполнены. Характер перемещений узлов также соответствует представлению о поведении ЗРС рамы под действием сил.

 


Понравилась статья? Добавь ее в закладку (CTRL+D) и не забудь поделиться с друзьями:  



double arrow
Сейчас читают про: