Пособие по МКЭ (864300), страница 8
Текст из файла (страница 8)
25,а имеют вид:N1 = (1 − ξ )(1 − η )(1 − ζ );N 2 = ξ (1 − η )(1 − ζ );N 3 = ξη (1 − ζ );N 4 = (1 − ξ )η (1 − ζ );N 5 = (1 − ξ )(1 − η )ζ ;N 6 = ξ (1 − η )ζ ;N 7 = ξηζ ;N8 = (1 − ξ )ηζ .Функции формы в местной системе координат для узлов двадцатиузлового элемента при нумерации рис. 25,в имеют вид:N1 = (ξ − 1)(2ξ − 1)(η − 1)(2η − 1)(1 − ζ )(1 − 2ζ );N 2 = ξ (2ξ − 1)(η − 1)(2η − 1)(1 − ζ )(1 − 2ζ );N 3 = ξ (2ξ − 1)η (2η − 1)(1 − ζ )(1 − 2ζ );N 4 = (ξ − 1)(2ξ − 1)η (2η − 1)(1 − ζ )(1 − 2ζ );N 5 = −(ξ − 1)(2ξ − 1)(η − 1)(2η − 1)ζ (1 − 2ζ );50Основы метода конечных элементовN 6 = −ξ (2ξ − 1)(η − 1)(2η − 1)ζ (1 − 2ζ );N 7 = −ξ (2ξ − 1)η (2η − 1)ζ (1 − 2ζ );N 8 = −(ξ − 1)(2ξ − 1)η (2η − 1)ζ (1 − 2ζ );N 9 = −4(ξ − 1)ξ (η − 1)(2η − 1)(1 − ζ )(1 − 2ζ );N10 = −4ξ (2ξ − 1)(η − 1)(1 − ζ )(1 − 2ζ );N11 = −4(ξ − 1)ξη (2η − 1)(1 − ζ )(1 − 2ζ );N12 = −4(ξ − 1)(2ξ − 1)(η − 1)η (1 − ζ )(1 − 2ζ );N13 = 4(ξ − 1)(2ξ − 1)(η − 1)(2η − 1)(1 − ζ )ζ ;N14 = 4ξ (2ξ − 1)(η − 1)(2η − 1)(1 − ζ )ζ ;N15 = 4ξ (2ξ − 1)η (2η − 1)(1 − ζ )ζ ;N16 = 4(ξ − 1)(2ξ − 1)η (2η − 1)(1 − ζ )ζ ;N17 = 4(ξ − 1)ξ (η − 1)(2η − 1)ζ (1 − 2ζ );N18 = 4ξ (2ξ − 1)(η − 1)ηζ (1 − 2ζ );N19 = 4(ξ − 1)ξη (2η − 1)ζ (1 − 2ζ );N 20 = 4(ξ − 1)(2ξ − 1)(η − 1)ηζ (1 − 2ζ ).Применение местной системы координат позволяет упроститьалгоритмизацию процесса вычисления матриц конечных элементов.Однако в этом случае для вычисления матриц конечных элементовтребуется применение численного интегрирования.Таким образом, чтобы определить значение матриц элемента,необходимо найти значения подынтегральных выражений для этихматриц в точках интегрирования и провести суммирование.После получения матриц жесткости элементов и векторов нагрузок элементов и формирования глобальной матрицы жесткости конструкции и глобального вектора нагрузок конструкции производятсяучет условий закрепления конструкции и решение линейного уравнения.
Тем самым определяются перемещения всех узлов под нагрузкой.Затем по величинам перемещений узлов находят поля деформаций инапряжений в элементах и, следовательно, во всей конструкции.Как уже говорилось выше, общепринятая формулировка МКЭпредполагает отыскание поля перемещений и тем самым связана сминимизацией потенциальной энергии системы при отыскании узловых значений вектора перемещений. После того, как перемещениябудут определены, производится вычисление компонент тензоровдеформаций и напряжений.51Основы метода конечных элементовМатрицы призматического элемента определяем численным интегрированием выражений для их компонент.
Для этого при произвольной форме элемента вводятся местные координаты, связанные сэлементом: ξ , η , ζ .Узлы нумеруются так, как показано на рис 25. При этом местные координаты связаны с декартовыми посредством функций формы. Местные координаты имеют следующие значения:- координата ξ изменяется от нуля на поверхности, заданной узлами 1, 4, 8, 5 до единицы на поверхности, заданной узлами 2, 3, 7, 6;- координата η изменяется от нуля на поверхности 1, 2, 6, 5 доединицы на поверхности 4, 3, 7, 8 координата ζ изменяется от 0 наповерхности 1,2,3,4 до 1 на поверхности 5, 6, 7, 8.Для численного интегрирования рационально применять квадратуру Гаусса-Лежандра, которая в сравнении с другими методамипри той же точности требует меньше точек интегрирования.Для восьмиузлового изометрического элемента необходимо подве точки в направлении каждой из местных координат.
Значение интеграла получается суммированием значений подынтегральной функции в точках интегрирования, умноженных на весовые коэффициенты:111i =n j =n k =n∫ ∫ ∫ f (ξ ,η , ζ )dξ dη dζ = ∑ ∑ ∑ H i H j H k f (ξi ,η j , ζ k ),000(1.70)i =1 j =1 k =1где n - число точек интегрирования по одному из направлений; H i ,H j , H k - весовые коэффициенты для квадратуры Гаусса-Лежандра;ξi , η j , ζ k - координаты точек интегрирования.При двух точках интегрирования в каждом направлении (дляматрицы жесткости восьмиузлового элемента) весовые коэффициенты H i , H j , H k равны и имеют значения 0,125, а координаты точекинтегрирования ξi , η j , ζ k принимают значения 0,5 − 3 / 6 и0,5 + 3 / 6 , т.е.
0,211325 и 0,788675.Таким образом, чтобы определить значение матриц такого элемента, необходимо найти значения подынтегральных выражений дляэтих матриц в точках интегрирования и провести суммирование. Рассмотрим выражения для матриц, входящих в подынтегральные зависимости. Матрица упругих характеристик в случае изотропного материала и трехосного напряженного состояния имеет следующий вид:52Основы метода конечных элементов⎡γ ϕ ϕ 0 0 0 ⎤⎢ϕ γ ϕ 0 0 0 ⎥⎢⎥⎢ϕ ϕ γ 0 0 0 ⎥[ D] = ⎢(1.71)⎥,000G00⎢⎥⎢0 0 0 0 G 0 ⎥⎢⎥⎣0 0 0 0 0 G⎦причем упругие характеристики определяются соотношениямиµEϕ=;(1 + µ )(1 − 2µ )(1 − µ ) Eγ=;(1 + µ )(1 − 2 µ )E,G=2(1 + µ )где E - модуль упругости первого рода для изотропного материала;µ - коэффициент Пуассона; G - модуль сдвига.Для удобства записи, обозначим через H iξ частную производ-ную функции формы N i узла по переменной ξ .
Таким образом:∂N∂N∂NH iξ = i ; H iη = i ; H iζ = i ;∂ξ∂η∂ζ∂N∂N∂N(1.72)H ix = i ; H iy = i ; H iz = i ;∂z∂x∂yТогда матрица градиентов запишется в следующем виде:00 . . .⎤⎡ H ix⎢ 0 H0 . . .⎥iy⎢⎥⎢ 00 H iz . . .⎥[ B] = ⎢(1.73)⎥,0...HHix⎢ iy⎥⎢ 0 H iz H iy . . .⎥⎢⎥0 H ix .
. .⎦⎣ H izгде i = 1, 2, 3, ..., n - индекс узла (для тетраэдрального элемента n = 4,для восьмигранного - n = 8).Матрица жёсткости конечного элемента вычисляется интегрированием по объёму элемента матричного произведения [ B]T [ D][ B] .53Основы метода конечных элементовМатрица, получаемая в результате такого перемножения, длявосьмиузлового конечного элемента имеет размер 24 × 24, для тетраэдрального – 12 × 12. Но эту матрицу можно представить в виде матрицы размером 8 × 8 (для тетраэдрального элемента - 4 × 4), каждыйэлемент которой является подматрицей размера 3 × 3. При этом каждая подматрица соответствует связи узла номер i и номер j.
Подматрица произведения имеет следующую структуру:ijijij ⎤⎡ k11k12k13⎢ ijijij ⎥(1.74)kkk⎢ 21 2223 ⎥ ,⎢ k ij k ij k ij ⎥33 ⎦⎣ 31 32гдеijk11= γH ix H jx + GH iy H jy + GH iz H jz ;ijk12= ϕH ix H jy + GH iy H jx ;ijk13= ϕH ix H jz + GH iz H jx ;ijk 21= ϕH iy H jx + GH ix H jy ;ijk 22= γH iy H jy + GH ix H jx + GH iz H jz ;ijk 23= ϕH iy H jz + GH iz H jy ;ijk31= ϕH iz H jx + GH ix H jz ;ijk32= ϕH iz H jy + GH iy H jz ;ijk33= γH iz H jz + GH ix H jx + GH iy H jy .Выражение (1.74) является универсальным и не зависит от количества узлов в конечных элементах, аппроксимирующих трехмерное тело, но зависит от вида матрицы упругих характеристик (1.71).При получении значения тройного интеграла в выражении (1.70)следует произвести замену переменных.
При этом необходимо выяснить соотношения, связывающие переменные x, y, z с локальнымипеременными ξ , η , ζ , т.е. найти функции видаx = q (ξ ,η , ζ ) ; y = g (ξ ,η , ζ ) ; z = h(ξ ,η , ζ ) .Функции q, g, h должны быть непрерывны вместе со своимипроизводными 1-го порядка и устанавливать взаимно однозначное и вобе стороны непрерывное соответствие между точками области интегрирования пространства x, y, z и точками области интегрирования54Основы метода конечных элементовпространства ξ , η , ζ . Также необходимо, чтобы функциональныйопределитель (определитель матрицы Якоби) этих функций сохранялв области интегрирования постоянный знак. Матрица Якоби в трёхмерном случае имеет следующий вид⎡ ∂x ∂y ∂z ⎤⎢ ∂ξ ∂ξ ∂ξ ⎥⎢⎥∂∂∂xyz⎥.[I ] = ⎢(1.75)⎢ ∂η ∂η ∂η ⎥⎢ ∂x ∂y ∂z ⎥⎢⎥⎣ ∂ζ ∂ζ ∂ζ ⎦При перечисленных выше условиях можно находить значениеинтеграла (1.70), используя местные координаты:∫∫∫ f ( x, y, z)dx dy dz = ∫∫∫ f [q(ξ ,η, ζ ), g (ξ ,η, ζ ), h(ξ ,η, ζ )] I dξ dη dζ . (1.76)VVПоскольку для восьмиузлового элемента преобразование координат записывается через функции формы согласно выражению(1.69), то получаем частные производные переменных x, y, z по переменным ξ , η , ζ , дифференцируя выражение (1.69) и применяя обозначения (1.72):∂x i =n= ∑ H iξ xi ;∂ξ i =1∂x i =n= ∑ H iη xi ;∂η i =1∂x i =n= ∑ H iζ xi ;∂ζ i =1∂y i =n= ∑ H iξ yi ;∂ξ i =1∂y i =n= ∑ H iη yi ;(1.77)∂η i =1∂y i =n= ∑ H iζ yi ;∂ζ i =1∂z i =n= ∑ H iξ z i ;∂ξ i =155Основы метода конечных элементов∂z i =n= ∑ H iη z i ;∂η i =1∂z i =n= ∑ H iζ z i ;∂ζ i =1Для изопараметрических элементов функции формы для интерполяции непрерывной величины и функции преобразования координат одинаковы.
Для субпараметрических элементов - различны, поскольку при определении непрерывной величины применяется больше узлов, чем при определении геометрии. Значение определителяматрицы Якоби I в точках интегрирования при численном интегрировании выражения (1.76) вычисляем через подстановку значенийточек интегрирования в выражение для матрицы Якоби (1.75).Производные функций формы в местных координатах получаемприменением операции дифференцирования к выражениям (1.64).Производные функций формы по глобальным координатам, входящие в выражение (1.73), определяются через производные функцийформы по местным координатам с использованием обозначений выражения (1.72) следующим образом:⎧ H iξ ⎫⎧ H ix ⎫⎪⎪⎪−1 ⎪(1.78)⎨ H iy ⎬ = [ I ] ⎨ H iη ⎬ .⎪H ⎪⎪H ⎪⎩ iz ⎭⎩ iζ ⎭В выражении (1.78) матрица, обратнаяляется соотношением:⎡ ∂y ∂z ∂z ∂y ∂z ∂x ∂x ∂z⎢ ∂η ∂ζ − ∂η ∂ζ ∂η ∂ζ − ∂η ∂ζ⎢II⎢∂x ∂z ∂z ∂x⎢ ∂z ∂y ∂y ∂z−−⎢ ∂ξ ∂ζ ∂ξ ∂ζ∂ξ ∂ζ ∂ξ ∂ζ−1[I ] = ⎢II⎢⎢ ∂y ∂z ∂z ∂y∂z ∂x ∂x ∂z−⎢ ∂ξ ∂η ∂ξ ∂η ∂ξ ∂η − ∂ξ ∂η⎢II⎢⎢⎣56матрице Якоби, опреде∂x ∂y ∂y ∂x ⎤−∂η ∂ζ ∂η ∂ζ ⎥⎥I⎥∂y ∂x ∂x ∂y ⎥−∂ξ ∂ζ ∂ξ ∂ζ ⎥⎥ .
(1.79)I⎥∂x ∂y ∂y ∂x ⎥−∂ξ ∂η ∂ξ ∂η ⎥⎥I⎥⎥⎦Основы метода конечных элементовМатрицу жёсткости элемента получаем численным интегрированием по объёму элемента матричного произведения [ B ]'[ D ][ B ] сприменением выражений (1.75) – (1.79).Аналогичным образом производится нахождение матриц нагрузки для элементов, составляющих глобальный вектор нагрузки.Представленная выше последовательность нахождения матрицыжёсткости элемента имеет то преимущество, что матрица определяется постепенно, элемент за элементом. Таким образом, отсутствует необходимость сохранения всей матрицы жёсткости элемента, что длявосьмиузлового элемента требует 576 чисел, а для двадцатиузлового 3600 чисел.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
