1625914689-e957c8b7a8e4003fe3539e4e0e465a65 (532400), страница 30
Текст из файла (страница 30)
е. на действиявнешних сил периодического вида) следующим образом. Пусть нас интересует вклад в решение (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 → ω, т.
е. тело испытывает явление резонансапри приближении частоты колебаний p внешней силы к частоте собственных колебанийω этого тела.1226.1.4.Глава 6. Основы численных методов решения задач деформирования тел . . .Прямое пошаговое интегрирование уравнений движенияРассмотрим решение системы (6.43) методом пошагового интегрирования. Сегментвремени [0, T ] разбиваем на подсегменты (шаги) с границами 0, ∆t, 2∆t, ..., T . Для начального момента времени t = 0 векторы 0 U, 0 U̇, 0 Ü считаются известными. Векторперемещений ансамбля узловых точек 0 U и вектор скоростей этих точек 0 U̇ определяются из начальных условий (6.44), а вектор их ускорений 0 Ü определяется из уравненийдвижения (6.43), отнесенных к моменту времени t = 0.
Требуется найти векторы перемещений, скоростей и ускорений t U, t U̇, t Ü для всех остальных моментов времени t. Припошаговом интегрировании уравнений движения (6.43) предполагается, что векторы U,U̇, Ü известны для моментов времени 0, ∆t, . . . , t − ∆t, t. Требуется найти эти векторыв момент времени t + ∆t, т. е. на следующем шаге интегрирования.6.1.4..1Пошаговое интегрирование уравнений движения по явной конечноразностной схемеРассмотрим метод центральных разностей, который принадлежит к явным схемаминтегрирования уравнений движения.
В этом методе ускорения и скорости представляются в виде центральных разностей со вторым порядком точности. Схема метода центральных разностей является трехслойной, т. е. в разностные уравнения входят величины длятрех последовательных шагов времени: t − ∆t, t и t + ∆t. Векторы скоростей t U̇ и ускорений t Ü узловых точек аппроксимируются выражениями со вторым порядком точности:tÜ =1 t−∆t(U − 2 t U + t+∆t U),∆t2tU̇ =1(−t−∆t U + t+∆t U).2∆t(6.73)Подставляя выражения для вектора t Ü из (6.73) в (6.43), получим систему линейныхалгебраических уравнений для определения вектора t+∆t U:()121t+∆ttMU = R − K − 2 M t U − 2 M t−∆t U.(6.74)2∆t∆t∆tОтдельно рассмотрим аппроксимации в момент времени t = 0. Решая систему (6.73) относительно векторов t−∆t U, t+∆t U, получим в «законтурной» точке, соответствующей моменту времени −∆t, следующее представление вектора перемещений U:∆t2 0Ü.(6.75)2Это выражение позволяет определить левую часть формулы (6.74) в момент времени t = 0и найти вектор перемещений ∆t U из этой формулы.Введем эффективную матрицу масс M̂ и эффективный вектор внешних сил R̂:()121ttM̂ ≡M, R̂ ≡ R − K − 2 M t U − 2 M t−∆t U.(6.76)2∆t∆t∆t−∆tU = 0 U − ∆t 0 U̇ +В этих обозначениях формула (6.74) переписывается следующим образом:M̂ t+∆t U = t R̂.(6.77)Примечание.
Если вместо прямого интегрирования рассматривается решение динамических уравнений методом разложения по формам собственных колебаний, то матрицамасс M ансамбля узловых точек заменяется на единичную матрицу I, матрица жесткостиTK – на матрицу Ω2 , а вектор внешних сил R – на вектор Φ̃ R (см. уравнения (6.43) и(6.64)).Алгоритм вычислений по явной схеме состоит из двух этапов. На первом этапепроводятся начальные вычисления:6.1. Применение метода конечных элементов к решению задач .
. .1231. Формирование матриц жесткости K и масс M.2. Введение векторов начальных значений 0 U, 0 U̇, 0 Ü.3. Задание шага ∆t < ∆tcr (∆tcr – размер шага, при превышении которого вычислительный процесс становится неустойчивым).4. Определение вектора−∆tU из формулы (6.75).5. Формирование эффективной матрицы масс M̂ с использованием формулы (6.76).6. Триангуляризация эффективной матрицы масс M̂ = LDLT (этот шаг не нужен вслучае модального интегрирования).На втором этапе проводятся вычисления в циклах, соответствующих шагам повремени. На каждом шаге выполняются три операции:1. Вычисление вектора t R̂ из формулы (6.76).2. Решение системы уравнений методом Краута (см. п. 6.1.2):(LDLT )t+∆t U = t R̂.(6.78)3.
Нахождение векторов t U̇, t Ü по формулам (6.73), затем переход к следующему шагу(t → t + ∆t).Этот алгоритм пошагового интегрирования уравнений движения (6.43) можно использовать в том случае, если использован вариант формирования совместной матрицымасс M (тогда эта матрица имеет тот же профиль, что и матрица жесткости K).Рассмотрим некоторые модификации этого алгоритма.Модификация 1.
Пусть матрица масс M – диагональная, что является допустимым упрощением при решении задач высокоскоростного деформирования тел с мелкойконечно-элементной сеткой. В этом случае система (6.77) распадается на NEQ независимых уравнений следующего вида:1mi (i = 1, NEQ),(6.79)∆t2где mi – значение i-го диагонального элемента матрицы M.
В этом случае отпадает необходимость в триангуляризации эффективной матрицы масс M̂, а система уравнений (6.78)заменяется решением NEQ уравнений (6.79).Модификация 2. Вместо того, чтобы сначала найти матрицу жесткости ансамбляузловых точек K, а потом определить произведение K t U в (6.76), можно при вычисленииэтого произведения в (6.76) использовать следующее его представление (см. (6.34)):m̂i t+∆t Ui = t R̂i ,tm̂i ≡K U=K∑K̃j t U.(6.80)j=1Такое представление произведения K t U позволяет определять только матрицы жесткостиэлементов K̃e и, таким образом, обойтись без хранения в памяти компьютера матрицыжесткости ансамбля узловых точек K.Модификация 3. Можно даже обойтись без вычисления матриц жесткостей элементов.
Справедливы следующие равенства:∫∫∫∫et eeTet eeTet eeTtK U = ( B CB dV ) U = B CB U dV = B C ε̄ dV = Be T t σ̄ dV.VeVeVeVe(6.81)124Глава 6. Основы численных методов решения задач деформирования тел . . .Введем вектор внутренних сил элемента:∫t eF ≡ Be T t σ̄ dV.(6.82)VeЭтот вектор вычисляется с использованием квадратурных формул Гаусса – Лежандра.Определим также вектор внутренних сил ансамбля узловых точек:tF≡K∑tjF̃ =j=1K∑Aj T t Fj .(6.83)j=1Тогда из (6.81)–(6.83) получим следующее равенство:ttK U= F=K∑j=1jF̃ =K∑Aj T t Fj .(6.84)j=1Таким образом, для вычисления произведения K t U в (6.76) достаточно вычислить вектор внутренних сил ансамбля узловых точек t F суммированием векторов внутренних силэлементов t F̃e .Достоинство явной схемы интегрирования уравнений движения (6.43) состоит в еепростоте и в экономном расходовании памяти компьютера. Однако у этой схемы есть инедостаток, который связан с ограничением шага интегрирования ∆t (т.
е. явная схема является условно устойчивой): вычисления устойчивы по отношению к ошибкам округленияв том случае, если выполнены неравенства:∆t ≤ ∆tcr = Tmin /π,(6.85)где Tmin – наименьший период собственных колебаний, т. е. Tmin = 2π/ωNEQ , гдеωNEQ – наибольшая собственная частота ансамбля узловых точек.
В том случае, если сеткаравномерна, критерий принимает более простой вид (критерий Куранта): ∆t < ∆l/c.Здесь ∆l – шаг сетки, c – скорость распространения возмущений в твердом теле (скорость√звука). Для одномерного случая (например, для стержня) скорость звука равна c = E/ρ.6.1.4..2Пошаговое интегрирование уравнений движения по неявной схемеНьюмаркаВ отличие от явной схемы, неявная схема Ньюмарка является двухслойной, так как вразностные уравнения входят величины не для трех последовательных моментов времени(как это было для явной центрально-разностной схемы), а для двух: t и t + ∆t. Принимаются следующие представления векторов перемещений и скоростей узловых точекансамбля в момент времени t + ∆t:t+∆tU̇ = t U̇ + [(1 − δ)t Ü + δ t+∆t Ü]∆t,t+∆tU = U + U̇∆t + [(1/2 − α) Ü + αttt(6.86)t+∆t2Ü]∆t .Здесь α и δ – постоянные параметры.Рассмотрим значения параметров α = 1/4, δ = 1/2.
В этом случае равенства (6.86)принимают следующий вид:1U̇ = t U̇ + (t Ü + t+∆t Ü)∆t,21t+∆tU = t U + t U̇∆t + (t Ü + t+∆t Ü)∆t2 .4t+∆t(6.87)6.1. Применение метода конечных элементов к решению задач . . .125Введем обозначение среднего ускорения:1a ≡ (t Ü + t+∆t Ü),2(6.88)используя которое, получимt+∆tU̇ = t U̇ + a∆t,t+∆tU = t U + t U̇∆t + a∆t2.2(6.89)Формулы (6.89) являются точными при постоянном значении вектора ускорений Ü напромежутке времени [t, t + ∆t].Уравнение (6.43), записанное в момент времени t + ∆t, принимает следующий вид:M t+∆t Ü + K t+∆t U = t+∆t R.(6.90)Введем следующие обозначения:a0 =1,α∆t2a1 =1,α∆ta2 =1− 1,2αa3 = ∆t(1 − δ),(6.91)a4 = δ∆t.Из второго соотношения (6.86) находим вектор t+∆t Ü и, считая его известным, из первогосоотношения (6.86) находим вектор t+∆t U̇.
В результате имеемt+∆tÜ = a0 (t+∆t U − t U) − a1 t U̇ − a2 t Ü,Подставляя выражение дляt+∆tt+∆tU̇ = t U̇ + a3 t Ü + a4 t+∆t Ü.(6.92)Ü из (6.92) в (6.90), получимK̂ t+∆t U = t+∆t R̂.(6.93)Здесь введены эффективные матрица жесткости K̂ и вектор внешних силK̂ ≡ K + a0 M,t+∆tt+∆tR̂:R̂ = t+∆t R + M(a0 t U + a1 t U̇ + a2 t Ü).(6.94)Алгоритм вычисления по схеме Ньюмарка состоит из двух этапов. На первом этапепроводятся начальные вычисления:1. Формирование матриц жесткости K и масс M.2.
Определение векторов начальных значений 0 U, 0 U̇, 0 Ü.3. Выбор шага ∆t, параметров α, δ. Если δ ≥ 0, 5 и α ≥ 0, 25(0, 5 + δ)2 , то схемаНьюмарка является безусловно устойчивой. При стандартных значениях δ = 0, 5, α =0, 25 имеем схему второго порядка точности. После выбора параметров производитсявычисление констант a0 , a1 ,.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
