Стержневой конечный элемент

[Avksent 10]

Здравствуйте!

 

Решил опробовать свои силы в МКЭ чуть глубже, а именно написать свою програмку для конечных элементов.

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

 

Однако пришла беда: матрица жесткости элемента получается странной, а именно - её определитель равен нулю

Что бы я не делал, но выявить причину мне не удалось 

 

Вот, обращаюсь к старейшинам за подсказкой, помощью и поддержкой.

 

Итак, собственно что я делал для вывода матрицы жесткости:

 

1. Рассмотрел задачу защемленного на конце (z=0) стержня, где z - осевая координата от 0 до L

*тут должна быть система уравнений теории стержней

 

2. Предположил, что объёмные нагрузки отсутствуют.

 

3. Собственно решил систему уравнений и получил решения типа:

 

u_i(L) = Сумма(Q_j*C_j)

 

где u_i(L) - соответствующее перемещение или поворот 

Q_j - Соответствующие силы или моменты

C_i - Какие-то константы, зависящие от свойств стержня

 

4. Очевидно, что всю эту красоту можно записать в матричном виде:

 

Uz = Sz Fz

 

Uz = [u₁(L),u₂(L),...u₆(L)]^T

 

Fz = [Q₁,Q₂,...M₃]

 

Sz = [...]_6х6- матрица размерностью 6х6

 

Стоит отметить, что часть элементов S равны нулю:

 

5. Обращая S (матрица податливости по-логике), получим матрицу жесткости для элемента с защемленным концом:

 

Kz = Sz^-1

 

6. Тогда потенциальная энергия элемента, обусловленная его деформацией, может быть записана в виде:

 

П^(e) = 0.5 Uz^T*Kz*Uz

 

7. Далее необходимо найти матрицу жесткости в общем виде, для чего полагаем все перемещения и повороты на конце z = 0 произвольными. Вспоминаем, что потенциальная энергия плевать хотела  инвариантна по отношению к жестким движениям, т.о. можно двигая элемент, как твёрдое тело, прийти к нулевым перемещениям и поворотам на конце z = 0, при этом значение П не изменится.

 

Тогда произвольным перемещениям и поворотам на концах

 

U₀ = [u₁(0),u₂(0),...u₆(0)]^T  и  U_L = [u₁(L),u₂(L),...u₆(L)]^T

 

Соответствует следующее значение Uz

Uz = U_L + Xz*U₀

где Xz = 

Или в другой записи через вектор узловых перемещений и поворотов элемента:

 

U^(e) = [u₀, U_L ]^T 

X = [Xz, Ez], Ez - единичный тензор 6х6

Uz = X*U^(e)

 

8. Вспоминаем о потенциальной энергии:

 

П^(e) = 0.5 Uz^T*Kz*Uz = 0.5 (X*U^(e))^T*Kz*(X*U^(e)) = 0.5 U^(e)T*K*U^(e)

 

9. Т.о. матрица жесткости элемента:

 

K = X^T*Kz*X

 

10. А её определитель, при определённых параметрах стержня равен нулю... какая жаль

 

 

 

Добавлю: определители Sz и Kz нулю не равны.

Изменено 20 февраля 2015 пользователем Avksent

[Борман 2 375]

Однако пришла беда: матрица жесткости элемента получается странной, а именно - её определитель равен нулю

Так оно и должно быть. Остальное почитаю...

[Avksent 10]

Так, судя по ответу Бормана мне пора немного отдохнуть и перечитать Зенкевича, Сегерлинда.

...

Так, не усну, пока не задам вопрос:

Пусть есть балка, разбитая на один КЭ, пусть она заделана в нуле... Ну конечно!!! Он же не закреплён никак у меня.

Спасибо!

 

Глупо создавать было тему...

Изменено 20 февраля 2015 пользователем Avksent

[Avksent 10]

Решил не создавать новую тему. Продолжу в этой ветке.

 

Итак, я дописал програмку, стал её проверять.

 

Первая проверка - статика. Искал решение уравнений вида: K.X = F, где K - матрица жесткости, X - столбец перемещений и поворотов, F - столбец узловых сил и моментов.

Различные решения сравнивал с аналитическими. Собственно, всё сошлось, в том числе и кручение стержня вокруг оси.

 

Затем решил найти собственные частоты. Соответственно имеем: A2.X'' + K.X = 0, A2 - матрица инерции. Решал в пакете WolframMathematica, поэтому код выглядит так (чтобы выдать частоту в Герцах):

Sqrt [ Eigenvalues [ { K, A2 } ] ] / ( 2*Pi )

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

 

Для примера: защемлённая балка (опоры не абсолютно жесткие). Длина: 1 м. Диаметр (сечение - круг): 0.02 м. Сталь. По длине 1 конечный элемент. (можно больше, но для примера одного мне кажется достаточным).

 

Результат собственных частот: {179663., 179663., 103382., 14616.5, 10284.4, 10284.4, 1391.05, 139.678, 139.678, 14.1885, 14.1885, 8.62927}

Если смотреть справа налево: вторая и третья СЧ - изгибные частоты (хорошо совпадают), шестая (1391.05) - продольные колебания (ошибка - 10%, при 4 КЭ ошибка 0.5%).

А вот дальше начинаются проблемы: первая частота (8.629) - имеет форму крутильных колебаний. Более того, при увеличении количества КЭ появляется много нереальных частот, имеющих крутильную форму.

Видимо ошибка содержится в матрице A2, её вид я выложу чуть позже. А пока без её вида, можете ли подсказать: куда копать?

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Отвечаю сам себе: проблема была с А2

невнимательность

Изменено 5 марта 2015 пользователем Avksent