Шабров Н.Н. - Метод конечных элементов в расчётах деталей тепловых двигателей (1061803), страница 11
Текст из файла (страница 11)
е. конечных элементов с линейной или простейшей из возможных аппроксимаций для искомых функций. Однако уже анализ решения задачи о растяжении стержня под действием собственного веса (см. рис. 1.5) показывает, что использование симплекс-элементов в этом случае не дает удовлетворительных результатов. При этом возникает естественное желание увеличить порядок интерполяционного поли- нома для перемещений. Например, от линейной аппроксимации перемещений в элементе перейти к квадратичной, чтобы при дифференцировании перемещений получить линейно изменяющиеся в элементе напряжения.
Можно повторить всю процедуру МКЭ для квадратичного интерполяциониого полинома в виде, аналогичном (1.2) (2О и(е) (х) = [1 х лз) и, (4.1) с(2 т. е. построить исходя из (4.1) функции формы для так называемого квадратичного конечного элемента и соответствующие матрицы и векторы. Но если такой путь и приемлем для решения одномерных задач, то для задач большей размерности ап является 69 менее удобным хотя бы из-за потери в общности алгоритма вычисления соответствующих интегралов, о чем упоминалось ранее. Введение локальной системы координат, связанной с конечным элементом, позволяет во многих отношениях стандартизировать вычислительные процедуры в общем алгоритме МКЭ.
Используя технику функций форм, определенных в локальной системе координат [см. ф-лу (2.54)), можно установить связь между глобальными и локальными координатами [8 ). Рис. 4Л. Линейный (а), квадратичный (б) иаонараметрические одномерные конечные элементы и субпараметрический элемент (в): СЧ вЂ” узлы, используемые для интерполяции перемещений; К вЂ” узлы для интерполяции формы элементе ой' а с) и' а бх д$ дх бЕ ' (4.6) 70 Вернемся к задаче о растяжении стержня под действием сил собственного веса и рассмотрим ее решение с новых позиций, используя один одномерный квадратичный конечный элемент в локальной системе координат (рис.
4.1, б) с постоянными пределами изменения независимой переменной — 1 ( $ ( 1. Такой элемент имеет три узла с одной степенью свободы в каждом из них. Функции формы для него в локальной системе координат имеют вид Ж; — — ф(1 — Е); У; = (1+~)(1 — $); Уй = ~ (1+Ц. (4.2) Приемы построения функций форм для различных конечных элементов описаны в работе [19). Интерполяционное соотношение для перемещений в элементе может быть представлено и " (Е) = [А' К) А(~ ($) А" (Б)) рм). (4.3) Аналогичным образом получим интерполяционное соотношение для глобальной координаты хСе) (Е) = К © Уу. (Е) Уй (ЕН Х('), (4.4) х; О Х(е) — х; = 1 (4.5) хй 21 Далее, следуя соотношению (1.18), вычислим производные от функций форм по глобальной координате. Поскольку теперь функции формы определены в локальной системе координат (4.21, можно написать Производная с[х/с1$ определяет матрицу Якоби 1, которая в случае одномерной задачи вырождается в число.
В соответствии с (4.4) можно написать 1, — [,'(( В/. 11/ ] 1~(с) ((Х д) (4.7) Решая (4.6) относительно (1Л/„/с)х с учетом (4,7), получим <Иа 1 — 1 с))( а ()х Нетрудно заметить, что матрица градиентов для рассматриваемой задачи будет иметь вид В(') = .1 — ' — [/)/1 Л'/ 1)1),). 11 (4.9) Определим матрицу Якоби для конкретного рассматриваемого конечного элемента.
Получим (4.8) О 3 — [ — — д- 1 — 2!в 21 (4.1О) Тогда (4. 11) В'~ = —, [ — — д-1 — 2! — -';1~ 1 — — +5 2 К(е] АŠ— — д-1 — 21 — д-2~ д! = 1 1 2 2 — 2$ 1 — +$ 2 — 8 1 16 — 8 — 8 7 7 АŠ— 8 61 1 (4. 13) Далее вектор узловых сил конечного элемента запишем 1 Г~(' = А ~ 1Ч" с( (1х = А ~ М" сд (1е1 Л (1$. 1(е) — 1 (4.14) 71 Матрица жесткости квадратичного конечного элемента согласно (1.24) будет представлена 1 К(') = А ) В(д)'ЕВ(') с1х = А ) В(')'ЕВ(2) с1е13 (11. (4.12) 1(д) — 1 В формуле (4.12) установлены пределы интегрирования и выполнено преобразование элемента длины в локальной системе координат.
Подставляя в (4.12) выражения для В(') из (4.11) и выполняя интегрирование, будем иметь 1 Подставляя в (4.14) выражения для Й(') из (4.2), после интегрирования получим ! Ра'» = А61 (4.15) — ! Поскольку рассматривается один конечный элемент, для формирования системы уравнений МКЭ матричное сложение типа (1.30) не потребуется, и система уравнений будет иметь вид К (е) (е) ~ (е) ил и 7 — 8 1 и; 2 — — 8 16 — 8 и; = —, 8 . (4.16) 1 — 8 7 и!) 2 В систему уравнений (4.16) следует внести граничное условие и! = О. Как известно, это равносильно вычеркиванию первого столбца и первой строки системы уравнений (4.16). В итоге получим следующую модифицированную систему: 16 — 8 и! А)) 8 Решение системы уравнений (4.17) дает точное значение перемещений в узловых точках элемента и совпадает с решением формулы (1.49). Однако в отличие от кусочно-линейного решения с использованием симплекс-элементов квадратичная интерполяция полученного решения в соответствии с (4.3) дает точное решение в любой точке конечного элемента.
Далее, вычисляя напряжения в конечном элементе, можно в соответствии с (1.21) записать О з а 2 Е ~ ч =- ЕВ~'~ч'~ = — [.— — -(-! — 2! — -!- (] = Я(1 — ~), (4,18) что соответствует точному решению для напряжений в любой точке конечного элемента. Действительно окончательный результат в (4,18) устанавливает линейный закон изменения напряжений, При этом на свободном конце стержня, т. е, при $ = 1, имеем точное значение (т =-- О, а на защемленном краю стержня, т, е. при 5 = — 1, — также точное значение а = 2Я.
Факт совпадения приближенного решения с точным для всех точек конечного элемента не должен вызывать удивления, так как 72 в рассматриваемом случае аппроксимация неизвестных квадратичным полиномом (4,3) соответствует точному закону изменения перемещений, Примечательно, что один квадратичный конечный элемент дает точное решение задачи, в то время как два симплекс- элемента дают неудовлетворительное приближенное решение.
Очевидно, что с изложенных выше позиций можно определить и одномерный симплекс-элемент (рис. 4.1, а), рассмотренный в гл. 1. Функции формы для такого элемента будут иметь вид (19] 2 ( ~)' 7 2 ( +~)' (4.19) Интерполяционные соотношения для перемещений и глобальной координаты х будут следующими: и~ ~ © = [У' Д) У~(1))~ч ~ . (4.20) (Ь) [М' Я) Л~т(Е)] Х~ ~ > (4.21) где [" ] (4.22) Для матрицы градиентов в данном случае можно написать В" =-- .1 ~~ [ж; М~), (4.23) где .1 = ~ [М; Л'~)Х~'~. (4.24) (4.25) 73 Нетрудно показать, что решение рассмотренной выше задачи о растяжении стержня, представленного двумя симплекс-элементами, определенными в локальной системе координат соотношениями (4.20) — (4.24), приводят к системе уравнений (1.45). При определении квадратичного и линейного конечного элемента в локальной системе координат для интерполяции искомых функций и установления связи между системами координат использовались одни и те же функции формы.
Это видно из соотношений (4.3), (4.4) для квадратичного элемента и (4.20), (4.21) для линейного элемента. Такие конечные элементы и рассмотренные ранее симплекс-элементы являются изопараметрическими элементами. Как указано в работе [19!, совпадение функций форм для интерполяции искомых функций и установления связи между координатами не является обязательным. В нашем примере интерполяционное соотношение (4.4) является квадратичным только формально, в то время как фактически оно является линейным, т.
е. 0 х<пщ [ — ~~ — ц ~~~ц(~ ц ~ ~~ ~ ~) ! -~~1)ц. 21 в вершинах четырехугольника, а последние четыре — на сере- динах сторон (рис. 4,2, б). Функции формы для такого элемента даются следующими соотношениями: Л(((5 Ч)-= — (1+ Б,)(1+ Ч(Ч) Я(+ ЧЧ; — 1)~ ((= 1, 2, 3, 4); (4.29) Л() Я, Ч) =- ~ (1 — ~')(1+ ЧЧ;) (! = 5, 7); 1 (4,3О) ~л„(е, Ч) = — (1 — Ч')(1 + Д(,) (А = — 6, 8). (4.31) Интерполяционное соотношение для перемещений в элементе будет иметь вид (4.32) (4.33) где х„и у„— глобальные координаты узлов элемента. Интерполяционному соотношению (4.33) можно придать другой вид: х(5 Ч) 1'(( ) (~ Ч) Х( ) (4.34) (4.35) Здесь И(') Д, Ч) — матрица-строка функций форм, а Х() и т'(е)— векторы узловых значений глобальных координат, т.
е. х1 х2 11(е) Х(е) (4.36) х!и 75 где т — число узлов конечного элемента (в данном случае равно 4 или 8); и и о — компоненты вектора перемещений в узловой точке с номером и. Аналогично строится интерполяционное соотношение для формы конечного элемента, которое устанавливает связь между глобальными декартовыми и локальными координатами По аналогии с (2,63) составим выражение для блока матрицы градиентов В (е) д — У дх 0 В~" = (4.37) Однако теперь функции формы в (4.37) определены в локальной системе координат в соответствии с соотношениями (4.29) — (4.31) или (4.28) и вычисление производных по глобальным координатам требует предварительных математических преобразований.
Эти преобразования аналогичны тем, которые были проделаны ранее в случае одномерного элемента, и необходимы для установления связи между производными функций форм по глобальным координатам и производным тех же функций по локальным координатам. Преобразование локальных производных в глобальные осуществляется соотношением д — Л' дх д д~ <е (4.38) д д ду д д дп где 1 — матрица Якоби, д — х д~ д — Д (4.39) д д д х д д дп дт~ Подставляя в (4.39) соотношения (4.34) и (4.35), получим д — Л~, д~ д д [ Х1 ) 71 ) 1 (4. 40) й~2 Лт д д 2 д т — ! — 1 76 Рассмотрим выражение для блока матрицы жесткости четырехугольного изопараметрического элемента, которое, очевидно, по форме совпадает с (2.67) К а) = ~ В~,') НВ~~'й) (а, ~1 = 1, 2, ..., т).
(4.41) А (е) Поскольку компоненты матрицы градиентов В~') определены в локальной системе координат при помощи соотношения (4.38), интеграл (4.41) приводится к удобной для вычисления форме: 1 1 К а = ~ ) В~ ) НВа) <1е1! Л !с5с1т1. (4.42) Подынтегральная функция в (4.42) в отличие от аналогичной функции в интеграле (2.67) для симплекс-элемента не является постоянной и зависит от локальных координат. Следует отметить, что в данном конкретном случае можно вычислить интеграл (4.42) в явном виде. Однако, как указывается в 1191, это скорее является исключением, чем правилом. Поэтому обычно используются ме- тоды численного интегрирования, среди которых наиболее рас- пространенным является метод Гаусса. Применяя одномерный метод численного интегрирования Га- усса по каждой координате, для интеграла (4,42) можно записать и л К!'в! = ~, ~ В„" НВв!'!1е1Л ~;.„, Н,Н;, ю=!(=! где $!, !~~ — точки интегрирования Гаусса по каждой координате; и — порядок интегрирования; Н; и Н, — весовые коэффициенты.