1626435462-e957c8b7a8e4003fe3539e4e0e465a65 (844208), страница 29

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 29)

Отметим, что принцип виртуальных работ можно также использовать для решения статических задач линейной теории упругости. Уравнение принципа виртуальныхработ записывается в виде (6.9). Выпишем выражения виртуальных работ внутренних ивнешних сил элемента:∫eδε̄T Cε̄ dV = δUe T Ke Ue , R̂Ve = δUe T ReV , R̂M= δUe T Me Üe , R̂pe = δUe T Rep .Ve(6.37)Здесь ввели симметричную матрицу масс элемента∫eM = HT Hρ dV.V(6.38)eИспользуя вектор неизвестных U ансамбля и выражения виртуальных работ (6.37), получим из (6.9) приближенное уравнение принципа виртуальных работ:TδUK∑(K̃ U + M̃ Ü −jjR̃jV−R̃jp )j=1−Kc∑Rik δuki = 0,i = 1, 2, 3.(6.39)k=1Вводим симметричную положительно определенную матрицу масс ансамбляM=K∑M̃j ,(6.40)j=1где матрица M̃e определяется аналогично матрице K̃e :δUe T Me Üe = δUT Ae T Me Ae Ü = δUT M̃e Ü,M̃e ≡ Ae T Me Ae .(6.41)С учетом введенных матриц K, M и вектора R аппроксимированное уравнение балансавиртуальных работ (6.39) запишем в следующем виде:δUT (KU + MÜ − R) = 0.(6.42)118Глава 6.

Основы численных методов решения задач деформирования тел . . .В силу произвольности вектора δU из (6.42) получим систему обыкновенных дифференциальных уравнений:MÜ + KU = R,(6.43)к которой добавляем начальные условия:U = U0 , U̇ = V0при t = 0.(6.44)Здесь U0 , V0 – заданные векторы (размерности NEQ) перемещений и их скоростей ансамбля узловых точек.6.1.2.Решение статических уравненийСистема линейных алгебраических уравнений (6.35) имеет матрицу K, обладающуюзамечательными свойствами. Во-первых, эта матрица симметрична. Во-вторых, она имеетпрофильный вид хранения элементов, т.

е. ненулевые элементы матрицы сосредоточеныоколо главной диагонали (см. (6.36)). В-третьих, эта матрица положительно определена всилу того, что справедливо неравенство:∫(6.45)ε̄T Cε̄ dV = UT KU ≥ 0.VПри подходящем задании граничных условий (6.7), когда перемещение тела как жесткогоцелого невозможно, равенство нулю в (6.45) достигается только в том случае, если U = 0,что и означает положительную определенность матрицы K.Отмеченные свойства симметричной матрицы K существенно облегчают решение системы (6.35) по сравнению с решением этой системы уравнений с матрицей, не обладающейуказанными свойствами. Например, при решении этой системы прямым методом Гауссане требуется делать выбор ведущего элемента.

Метод Гаусса решения системы уравнений(6.35) предполагает разложение матрицы K в виде K = LQ, где L – нижняя треугольнаяматрица, а Q – верхняя треугольная матрица. Для симметричной матрицы существует вариант метода Гаусса – метод Краута, когда производится разложение матрицы Kв виде K = LDLT . Здесь D – диагональная матрица с положительными компонентамиd11 , . . .

, dNEQ NEQ , а L – нижняя треугольная матрица, имеющая cтруктуру следующеговида:1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0  1 0 0 0 0 0 0 0  0 0 0 1 0 0 0 0 0 0  0 1 0 0 0 0 0 .L=(6.46)000010000 0 0 0 0 0 1 0 0 0  0 0 1 0 0  0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1Затем решается система:(LDLT )U = R⇔L(DLT U) = R⇔LV = R.(6.47)Решение последней системы уравнений в (6.47) соответствует выполнению прямого ходаметода Гаусса – нахождению вектора V. Затем решается система:DLT U = V⇔LT U = D−1 V.(6.48)6.1.

Применение метода конечных элементов к решению задач . . .119Решение последней системы уравнений в (6.48) соответствует выполнению обратного ходаметода Гаусса – нахождению вектора U, т. е. определению вектора перемещений ансамбля узловых точек. После этого в каждом элементе с помощью формулы (6.29) находитсявектор перемещений элемента Ue . Далее, по формулам (6.20), находятся компоненты тензора деформаций в элементе.

Компоненты тензора напряжений находятся из соотношений(6.3).Матрицы и векторы элементов находятся численным интегрированием выражений,представленных в п. 6.1.1. Это численное интегрирование проводится с использованиемформул Гаусса – Лежандра. Из теории МКЭ известно, что компоненты вектора перемещений находятся наиболее точно в узловых точках ансамбля, а компоненты тензорадеформаций, а следовательно, и тензора напряжений – в точках интегрирования внутриэлементов.Еще одним известным методом решения системы (6.35) является метод сопряженных градиентов, который принадлежит к итерационным методам решения этой системыалгебраических уравнений.

Итерационные методы решения таких систем уравнений являются альтернативными методами по отношению к прямым методам (например, методуГаусса) решения этих систем уравнений. При EBE (element by element) реализации метода сопряженных градиентов матрица K ансамбля узловых точек не хранится в памятикомпьютера, что позволяет решать задачи большой размерности (∼ 106 уравнений) наперсональных компьютерах.6.1.3.Решение динамических уравнений методом разложения решения по формам собственных колебанийРассмотрим уравнения собственных колебаний тела (т.

е. колебаний при отсутствиивнешних сил):MÜ + KU = 0,(6.49)которые получаются из уравнений (6.43) приравниванием вектора внешних сил нулевомувектору. Представим решение уравнения (6.49) в следующем виде:U = Φ sin ω(t − t0 ).(6.50)После подстановки этого выражения для вектора U в (6.49) получим обобщенную задачупо определению собственных значений/векторов:(K − ω 2 M)Φ = 0.(6.51)Решением полученной системы (6.51) являются собственные пары2(ω12 , Φ1 ), (ω22 , Φ2 ), . . . , (ωNEQ, ΦNEQ ),(6.52)где ωi2 – собственные значения, а Φi ̸= 0 (i = 1, NEQ) – М-ортогональные собственныевекторы, для которых выполнены следующие равенства:ΦTi MΦj = δij (i, j = 1, NEQ).(6.53)Величины ωi и Φi имеют механический смысл частот и форм собственных колебаний.Собственные пары в (6.52) ранжируются таким образом, что для частот собственных колебаний выполнены следующие неравенства:0 ≤ ω1 ≤ ω2 ≤ ...

≤ ωNEQ .(6.54)120Глава 6. Основы численных методов решения задач деформирования тел . . .Если 0 = ω1 = ω2 = ... = ωI (1 ≤ I ≤ NEQ), тогда из (6.51) и того, что Φi ̸= 0, следует, чтоматрица K вырождена, т. е. det K = 0. Это означает, что тело не закреплено как жесткоецелое, а собственные векторы Φi i = 1, I соответствуют таким векторам перемещений U,при которых происходит движение тела без деформаций.

В этом случае или не существует решения статических уравнений (6.35), или существует их неединственное решение,однако решение динамических уравнений можно найти по стандартным алгоритмам, приведенным далее.Задачу (6.43), (6.44) можно решить прямым интегрированием (см. п. 6.1.4). В этомслучае решение находится в арифметическом пространстве RNEQ c ортонормированнымбазисом:   100 1  0  0       e1 =  ...  , e2 =  ...  , .

. . , eNEQ =  ...  .(6.55)    0  0  0 001Решение U представляется в следующем виде:∑NEQU=Ui ei .(6.56)i=1Альтернативный способ решения задачи (6.43), (6.44) состоит в использовании метода разложения решения по формам собственных колебаний. Суть этого метода состоит в том, чтобы при нахождении решения в арифметическом пространстве RNEQ вместоортонормированного базиса e1 , e2 , . . . eNEQ использовать базис из М-ортогональных собственных векторов Φ1 , Φ2 , .

. . ΦNEQ .Введем матрицыω12 0 . . .0 0 ω2 . . .0 22Φ̃ ≡ [Φ1 , Φ2 , . . . , ΦNEQ ], Ω =  ..(6.57)......  . .... 20 0 . . . ωNEQИз (6.51), (6.57) получимKΦ̃ = MΦ̃Ω2 .(6.58)Из (6.53) имеем равенствоTΦ̃ MΦ̃ = I,(6.59)а из (6.58), (6.59) – равенствоTΦ̃ KΦ̃ = Ω2 .(6.60)Рассмотрим представление неизвестного вектора перемещений U в базисе М-ортогональныхсобственных векторов Φi (i = 1, NEQ):∑NEQU=Φi xi (t) = Φ̃X(t),(6.61)i=1где X – вектор коэффициентов разложения (6.61) видаXT = [x1 , x2 , ..., xNEQ ].(6.62)6.1. Применение метода конечных элементов к решению задач .

. .121Подставив представление решения (6.61) в уравнение (6.43), получим(6.63)MΦ̃Ẍ + KΦ̃X = R.TУмножая выражения в левой и правой частях равенства (6.63) на матрицу Φ̃ , получим,используя равенства (6.59), (6.60), систему обыкновенных дифференциальных уравненийотносительно величин xi (i = 1, NEQ) видаTẌ + Ω2 X = Φ̃ R.(6.64)Перепишем начальные условия (6.44), используя представление решения (6.61):Φ̃X(0) = U0 ,(6.65)Φ̃Ẍ(0) = V0 .TУмножая левые и правые части равенств (6.65) на матрицу Φ̃ M, получим, с учетомравенства (6.59), начальные условия для вектора X(t):TX(0) = Φ̃ MU0 ,TẌ(0) = Φ̃ MV0 .(6.66)Система (6.64) является несвязанной системой уравнений.

Она распадается на NEQ обыкновенных дифференциальных уравнений второго порядка:ri ≡ ΦTi R(t)(6.67)ẋi |t=0 = ΦTi MV0 ≡ vi0 .(6.68)ẍi (t) + ωi2 xi (t) = ri (t),с начальными условиями:xi |t=0 = ΦTi MU0 ≡ x0i ,Решение каждого уравнения можно получить либо численно, либо с помощью интеграла Дюамеля (для тела, в котором невозможны движения как жесткого целого, т. е. приусловии ωi > 0 (i = 1, NEQ)):1xi (t) =ωi∫tri (τ ) sin ωi (t − τ )dτ + αi sin ωi t + βi cos ωi t,(6.69)0где αi и βi определяются из начальных условий (6.68).Можно ограничиться нахождением приближенного решения U с учетом вклада врешение только нескольких нижних частот (пусть их число равно L) и заменить суммуот 1 до NEQ в (6.61) суммойL∑U≈Φi xi (t).(6.70)i=1Можно исследовать поведение системы на частотные воздействия (т.

е. на действиявнешних сил периодического вида) следующим образом. Пусть нас интересует вклад в решение (6.61) частоты i. Для краткости изложения опускаем нижний индекс i. Рассмотримрешение задачи:(6.71)ẍ(t) + ω 2 x(t) = R sin pt, x|t=0 = 0, ẋ|t=0 = 1.Из (6.69) и начальных условий в (6.71) получим()R/ω 21Rp/ω 3x(t) =sin pt +−sin ωt.1 − p2 /ω 2ω 1 − p2 /ω 2(6.72)Из (6.72) следует, что x(t) → ∞ при p → ω, т.

Характеристики

Список файлов книги