Переход к глобальным координатам

Матрица жесткости элемента записана в локальной системе координат (х, у, z).

Преобразование к глобальным координатам (х^’, у^’, z^’) необходимо для составления ансамбля элементов и затем соответствующих уравнений равновесия (рис.41).

EMBED Word.Picture.8

Рис.41

Узловые силы и перемещения преобразуются из глобальной в локальную систему координат с помощью матрицы [L]

где и

- косинус угла между осями х и х^’.

Для полной системы узловых сил элемента:

После преобразования матрица жесткости элемента в глобальных координатах принимает вид:

где

диагональная матрица, составлена из нескольких матриц [L], количество которых равно числу узлов элемента.

EMBED Word.Picture.8

Рис.42

Одну из осей локальной системы координат (х) направляем параллельно| глобальной оси х^’(рис.42). Для элемента i j k m определяем все направляющие косинусы. Направляющие косинусы оси х:

Направляющие косинусы оси у:

Направляющие косинусы для оси z:

Приведение заданной нагрузки к эквивалентным внешним узловым силам.

На элемент S действует нагрузка P_z(x, y), работа которой на виртуальных перемещениях, соответствует аппроксимирующей функции (28) определяется выражением:

(45)

Обобщенные силы определяются зависимостью:

(46)

Обобщенные силы являются эквивалентными внешними узловыми силами элемента S.

Нумерация и положительные направления для обобщенных сил и обобщенных координат q_j совпадают.

Вектор: является грузовым вектором элемента S.

Узел i (i = 1, 2,…N) в составе всей системы испытывает совместное действие внешних узловых сил, эквивалентных нагрузкам на объединенные элементы, содержащие узел i. Суммируя все обобщенные силы по каждому из направлений узла i (по 3 пластинок) получаем вектор внешней нагрузки узла i:

Блочный вектор компонентами которого являются векторы внешней нагрузки всех узлов системы, называется грузовым вектором системы.

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

где - силы со стороны элемента на узел.

Так как = [r] {q}, получим:

(47)

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

Чтобы исключить нулевые вариации, введем диагональную матрицу с N компонентами, записанными в том же порядке, что компоненты вектора перемещений. Номерам заданных перемещений соответствуют нулевые значения в матрице , остальные компоненты 1. Тогда с учетом правил транспортирования произведения матриц получим разрешающее уравнение равновесия метода конечных элементов:

.

Моделирование оболочечных конструкций

При моделировании методом конечных элементов (МКЭ) [1] широкого класса оболочечных систем, состоящих из монолитно соединенных подконструкций, имеющих срединные поверхности вращения, а также тонкостенных пространственных конструкций, криволинейных в плане, возникает проблема выбора формы конечных элементов. Использование наиболее простого плоского треугольного элемента не обеспечивает требуемого совпадения сетки триангуляции с линиями главных кривизн срединной поверхности указанных объектов. Поэтому представляется рациональным использование пластинчатых конечных элементов в форме равнобочной трапеции с узлами в углах. Принципиальные преимущества такого элемента при аппроксимации оболочек со срединными поверхностями вращения определяются тем, что границы этого элемента совпадают с линиями главных кривизн срединных поверхностей, что повышает точность итоговых результатов как в статических, так и в динамических задачах. Такие элементы естественно вписываются в сетку меридианов и параллелей оболочек со срединными поверхностями вращения и, получающаяся при измельчении сетки, многогранная поверхность достаточно полно аппроксимирует срединную криволинейную поверхность оболочки. Рассматриваемые конечные элементы испытывают суперпозицию изгиба и мембранного напряженного состояния.

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

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

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

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

Рис. 43. Мембранные и изгибные степени свободы элемента

Функции формы трапециевидного элемента в мембранном состоянии:

Функции формы, соответствующие изгибным степеням свободы узла i:

Матрица инертности пластинчатого элемента строится на основе общего выражения [1] для ее компонентов

где: e – порядковый номер элемента в ансамбле; v – объем элемента;

 – плотность материала.

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

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

Аналитические выражения для коэффициентов матрицы инертности в мембранном состоянии узла i имеют вид (здесь обозначено k = (b – a)/l):

Выражения для первых трех коэффициентов матрицы инертности изгибного состояния

( ):

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

Исследование частот и форм собственных колебаний упругих систем осуществлено на основе матричного уравнения , соответствующего обобщенной проблеме собственных значений. Метод конечных элементов, примененный к задачам такого типа, приводит согласно [7] к уравнению

,

где: [K] – матрица жесткости системы; [M] – матрица инертности;

– собственные значения.

Определение частот и форм выполнено в рамках неполной проблемы собственных значений [7]. Для ее решения выбран наиболее простой метод обратной итерации.

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

1. Собственные частоты колебаний прямоугольной упругой пластинки, защемленной по одной короткой стороне и свободной по трем другим (рис. 44).

Рис. 44. Сетка элементов защемленной прямоугольной пластинки: (a=2.5410^-2 м, b=2a, t=0.2510^-2 м, E=2.0210¹¹ H/м²,  = 0.3,  = 7798 кг/м³)

В таблице 1 приведены значения первых трех частот, найденные различными способами. При исследованиях использованы результаты решения рассматриваемой проблемы в среде программного продукта Pro/MECHANICA [8].

Таблица 1

Низшие частоты собственных колебаний прямоугольной пластинки

i, Гц	Аналитическое решение [9]	МКЭ, 4 треуг. элемента	МКЭ, 2 трапец. элемента	Pro/MECH.
	815	805	815	813
	3534	3632	3558	3373
	5121	5025	5121	5012

2. Собственные частоты колебаний косоугольной упругой пластинки (рис. 45), защемленной вдоль двух смежных сторон.

Рис. 45. Конечноэлементная сетка косоугольной пластинки: (a=610^-2 м, t=0.2510^-2 м, =26, E=2.0210¹¹ H/м²,  = 0.3,  = 7798 кг/м³)

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

Таблица 2

Низшие частоты собственных колебаний косоугольной пластинки

Трапециевидный элемент	1004	3641	4554	7155	10432	12069
Pro/ MECHANICA	1024	3820	4190	7467	10013	11196
Погрешность, %	2	5	8	4	4	8

3. Найдены собственные частоты и формы кольцевой пластинки (рис. 46), защемленной вдоль стороны большего радиуса.

Рис. 46. Конечноэлементная сетка защемленной кольцевой пластинки

Низшие частоты собственных колебаний, вычисленные аналитически [9] и полученные при конечноэлементном моделировании представлены в таблице 3.

Таблица 3

Низшие частоты собственных колебаний защемленной кольцевой пластинки

i, Гц	Аналитическое решение	Трапециевидный элемент	Треугольный элемент	Pro/ MECHANICA
	8126	7998	8479	8011
	14071	12644	13908	12453
	—	13053	14402	12453
	—	20922	21619	20630

Собственные формы, соответствующие найденным частотам показаны на рис. 47–50.

Рис. 47. Низшая собственная форма защемленной кольцевой пластинки для ₁

Рис. 48. Низшая собственная форма защемленной кольцевой пластинки для ₂

Рис. 49. Низшая собственная форма защемленной кольцевой пластинки для ₃

Рис. 50. Низшая собственная форма защемленной кольцевой пластинки для ₄

4. Тестирование матрицы инертности трапециевидного элемента проведено также при определении собственных частот цилиндрической оболочки, защемленной с двух сторон (рис.  51). Пример взят из [10], размеры указаны на рисунке.

Рис. 51. Защемленная цилиндрическая оболочка (размеры в м,  =1 кг/м³, E = 10⁴ Па,  = 0.2)

Решение сравнивалось с результатами исследований при помощи комплекса Pro-MECHANICA, а также конечноэлементного моделирования пластинчатыми элементами треугольной формы и элементом CAU2W30 (Табл. 4).

Таблица 4

Минимальная частота собственных колебаний защемленной цилиндрической оболочки

Низшая частота	Трапециевидный элемент	Треугольный элемент	Элемент CAU2W30	Pro/ MECHANICA
i, Гц	3.36	3.61	3.50	3.47

Расхождение результатов при определении низшей собственной частоты не превышает 5%.

5. Поведение цилиндрической оболочки с вырезом при собственных колебаниях рассмотрено на примере, описанном в [11]. На рис. 52 изображена исследуемая оболочка и приведены ее параметры.

Рис. 52 Цилиндрическая оболочка с вырезом (размеры в м,  =7800 кг/м³, E = 2.110¹¹ Па,  = 0.3)

Собственные частоты, полученные при конечноэлементном моделировании с использованием тестируемого элемента отличаются от соответствующих значений представленных в [11] не более, чем на 5%. Некоторые формы собственных колебаний показаны на рис. 53–56.

6. В качестве примеров использования трапециевидного конечного элемента рассмотрены задачи определения собственных частот и форм конической оболочки и оболочки в виде сектора тора. Оболочка толщиной 0.110^–2 м в виде прямого усеченного конуса длиной 1510^–2 м защемлена по стороне, имеющий радиус. 710^–2 м; сторона конуса меньшего радиуса 410^–2 м – свободна. В таблице 5 приведены собственные частоты, определенные численно разными авторами.

Рис. 53. Низшая собственная форма цилиндрической оболочки с вырезом для ₁

Рис. 54. Низшая собственная форма цилиндрической оболочки с вырезом для ₂

Рис. 55. Низшая собственная форма цилиндрической оболочки с вырезом для ₃

Рис. 56. Низшая собственная форма цилиндрической оболочки с вырезом для ₄

Таблица 5

Низшие частоты собственных колебаний защемленной конической оболочки

Трапециевидный элемент	1337	1337	1645	1645	1774	1774
Pro/MECHANICA	1324	1324	1786	1786	1825	1826
Погрешность, %	1	1	8	8	3	3

На рис. 57 показаны некоторые характерные формы собственных колебаний низших частот.

Оболочка в виде сектора кругового тора с углом раствора 45 имеет радиус поперечного сечения 710^–2 м, радиус образующей 2110^–2 м и толщину 0.110^–2 м. Материал оболочки имеет следующие характеристики: E = 2.110¹¹ Н/м²,  = 0.3,  = 7.8510³ кг/м³. Одно крайнее поперечное сечение защемлено, другое – свободно. В таблице 6 представлены результаты численного анализа.

Таблица 6

Низшие частоты собственных колебаний оболочки в виде сектора кругового тора, защемленной по одному крайнему поперечному сечению

i, Гц	Трапециевидный элемент	Pro/ MECHANICA	Погрешность, %
	569	558	2
	569	566	0.5
	952	1050	9

На рис. 58 изображены формы собственных колебаний оболочки.

^₁ ^₃ ^₅

Рис. 57. Низшие собственные формы конической оболочки

^₁ ^₃ ^₅

Рис. 58. Низшие собственные формы сектора торообразной оболочки

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

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

• аналитические выражения для компонентов матриц инертности повышают точность и быстроту их вычисления, по сравнению с численным интегрированием;

• обеспечивает возможность моделирования оболочек с произвольным изменением толщин элементов вдоль меридианов и параллелей срединной поверхности вращения;

• алгоритм генерации конечноэлементной сетки для оболочек со срединными поверхностями вращения легко поддается автоматизации;

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