1625914689-e957c8b7a8e4003fe3539e4e0e465a65 (532400), страница 32
Текст из файла (страница 32)
Схема применения модифицированного метода Ньютона – Рафсона показана на рис. 6.3, б. Недостатком этого метода по сравнению со стандартным являетсяболее медленная сходимость, к тому же итерационный процесс для этого метода чащеоказывается расходящимся по сравнению со стандартным методом.Cтандартный метод Ньютона – Рафсона является мощным инструментом решениязадач квазистатического деформирования, но он трудоемок в применении.
Модифицированный метод существенно снижает число операций на каждой итерации, но при этомухудшаются свойства сходимости. Кроме этих двух методов имеется семейство квазиньютоновых методов, которые по характеристикам сходимости и числу операций на однойитерации занимают промежуточное положение между стандартным и модифицированным методами Ньютона – Рафсона. Рассмотрим один из них – BFGS (Broyden – Fletcher –Goldfarb – Shanno) – квазиньютонов метод . Определим приращение вектора перемещенийδ (i) ≡ t+∆t U(i) − t+∆t U(i−1)(6.125)и приращение вектора несбалансированной нагрузкиγ (i) ≡ (t+∆t R − t+∆t F(i−1) ) − (t+∆t R − t+∆t F(i) ).(6.126)Квазиньютоновы уравнения записываются в следующем виде:t+∆tK̃(i) δ (i) = γ (i) .(6.127)Для BFGS (квазиньютонова) метода матрица (t+∆t K̃−1 )(i) в (6.127) находится по формуле:(t+∆t K̃−1 )(i) = A(i) T (t+∆t K̃−1 )(i−1) A(i) ,(6.128)где матрица A(i) размера NEQ × NEQ имеет следующий вид:A(i) = I + v(i) w(i)T ,(6.129)6.2.
Применение метода конечных элементов к решению задач теории . . .Rt+Dtt+DtR-t+Dt(1)t+DtFt+DtR-131(2)FRgtgKdgtdFt+Dtd(2)(3)(3)(2)(1)t+Dt(2)K(1)(1)K0t0Ut+Dtt+Dt(1)Ut+Dt(2)UUUРис. 6.4. Схема работы BFGS (квазиньютонова) методагде вводятся вектор-столбцы w(i) , v(i) вида√w(i)δ (i)≡ (i)T,δ γ (i)v(i) ≡ −t+∆t K̃(i−1) δ (i)δ (i)Tδ (i)T γ (i)− γ (i) .t+∆t K̃(i−1) δ (i)(6.130)В практических расчетах матрица (t+∆t K̃−1 )(i) в явном виде не вычисляется, вместо этогоопределяется вектор∆U(i) = (I + w(i−1) v(i−1) T ) .
. . (I + w(1) v(1) T ) t K−1 (I + v(1) w(1) T ) . . .. . . (I + v(i−1) w(i−1) T )(t+∆t R − t+∆t F(i−1) ). (6.131)Схема работы BFGS (квазиньютонова) метода приведена на рис. 6.4.Представленные выше итерационные методы можно ускорить при помощи процедуры линейного поиска. Для этого во всех методах на каждой итерации вместо определениявектора ∆U(i) находится вектор ∆Ū (вектор ∆Ū не меняется при выполнении линейного поиска), пробные значения вектора t+∆t U(i) , обозначаемые как t+∆t Ũ(i) , находятся поформуле:t+∆t (i)Ũ = t+∆t U(i−1) + β∆Ū.(6.132)Если в этом выражении положить β = 1, придем к обычным процедурам любого из используемых методов. Параметр β варьируется до тех пор, пока вектор невязки в направлении∆Ū не станет близок к нулю, т.
е. пока не будет выполнено следующее условие:∆ŪT (t+∆t R − t+∆t F̃(i) ) ∼ 0,где вектор t+∆t F̃(i) пересчитывается через вектор(6.133) достигается при выполнении условия:t+∆t(6.133)Ũ(i) . Приближенное равенство в∆ŪT (t+∆t R − t+∆t F̃(i) ) ≤ εs ∆ŪT (t+∆t R − t+∆t F(i−1) ),(6.134)где εs – заданный параметр. Считается, что оптимальным для параметра εs является значение εs = 0,5. Когда условие (6.134) оказывается выполнено, найденный вектор t+∆t Ũ(i)принимаем за вектор t+∆t U(i) .132Глава 6. Основы численных методов решения задач деформирования тел .
. .Для окончания итерационного процесса требуется задать критерии, по которым происходит переход на следующий шаг интегрирования уравнений квазистатического движения. Введем обозначение:U(i) ≡ t+∆t U(i) − t U.(6.135)Из (6.123) и (6.135) получим∆U(i) = U(i) − U(i−1) .Вводим евклидову норму вектора a = [a1 , a2 , . . .
, aN ]T следующим образом:√∥a∥2 ≡ a21 + a22 + . . . + aN1 .(6.136)(6.137)Контроль сходимости итерационных процессов можно осуществлять по трем критериям:1. Контроль сходимости по перемещениям∆U(i) 2≤ εD .(k)max ∥U ∥2(6.138)0≤k≤i2. Контроль сходимости по вектору невязки (вектору несбалансированной силы)t+∆tR − t+∆t F(i−1) 2≤ εF .(6.139)max∥τ R − τ −∆t F∥2τ =∆t,...,t,t+∆t3. Контроль сходимости по энергии∆U(i) T (t+∆t R − t+∆t F(i−1) )≤ εE .U(1) T (t+∆t R − t+∆t F)(6.140)В неравенствах (6.138)–(6.140) введены заданные положительные параметры εD , εF , εE .Обычно задаются значения εD = εF = 0,01, εE = 0,0001.Рассмотрим алгоритм решения квазистатических задач теории пластичности (дляBFGS (квазиньютонова) метода добавляем линейный поиск).Первый этап – начальные вычисления – состоит из двух действий:1.
Задание погрешностей вычисления εD , εF , εE и максимального числа шагов интегрирования kmax .2. Присвоение начального значения вектору 0 U.Второй этап – пошаговое интегрирование – состоит из следующих шагов (интегрирование идет по циклу, вначале полагаем k = 1):1. Касательная матрица жесткости t K вычисляется и разлагается на множители t K =LDLT .2. Вычисляется вектор невязки (эффективный вектор нагрузки, вектор несбалансированных внешних и внутренних сил) t+∆t R̂ = t+∆t R − t F.3. Решается система линейных алгебраических уравнений:(LDLT )U = t+∆t R̂,(6.141)в результате решения находится вектор приращений перемещений ∆U = t+∆t U(1) − t U.4. Проводим итерации равновесия, начиная с номера итерации i = 1, полагая,что U(1) = ∆U (здесь и далее знак «:=» обозначает присваивание значения некоторойвеличине).6.2.
Применение метода конечных элементов к решению задач теории . . .1334.1. i := i + 1.4.2. Находим (i − 1)-ю аппроксимацию вектора перемещений ансамбля узловых точек:t+∆tU(i−1) = U(i−1) + t U.(6.142)4.3. Находим (i−1)-ю аппроксимацию вектора несбалансированных внешних и внутреннихсил:t+∆t (i−1)R̂≡ t+∆t R − t F(i−1) .(6.143)4.4. Вариант для метода Ньютона – Рафсона. Если используем модифицированный методНьютона – Рафсона, то переходим к пункту 4.4.3.4.4.1. Находим матрицу касательной жесткости на i − 1-й итерации t+∆t K(i−1) .4.4.2. Разлагаем матрицу t+∆t K(i−1) на множители: t+∆t K(i−1) = LDLT .4.4.3. Решаем систему линейных алгебраических уравнений:(LDLT )∆U(i) = t+∆t R̂(i−1) ,(6.144)для модифицированного метода Ньютона – Рафсона используем разложение LDLTиз п.
1 второго этапа.4.4. Вариант для квазиньютонова метода с выполнением линейного поиска.4.4.1. Находим пробную коррекцию вектора приращения перемещений:∆Ū = (t+∆t K̃−1 )(i−1)t+∆t R̂(i−1) .(6.145)4.4.2. Осуществляем линейный поиск в направлении вектора ∆Ū, находим параметр β.4.4.3. Находим i-ю аппроксимацию вектора перемещений:t+∆tU(i) = t+∆t U(i−1) + β∆Ū.(6.146)4.4.4. Находим векторы:δ (i) = t+∆t U(i) − t+∆t U(i−1) ,γ (i) = t+∆t R̂(i−1) − t+∆t R̂(i) .(6.147)4.4.5.
С помощью векторов δ (i) , γ (i) находим новую аппроксимацию матрицы (t+∆t K̃−1 )(i) .4.5. Находим новое приращение вектора перемещений:U(i) = U(i−1) + ∆U(i) .(6.148)4.6. Делаем проверку на сходимость по критериям (6.138)–(6.140) (не обязательно по всем,можно делать проверку только по некоторым из них). Если критерии выполнены, то присваиваем ∆U := U(i) и переходим к п.
5. Если критерии не выполнены, то возвращаемсяк п. 4.1.5. Находим вектор перемещений в момент времени:t+∆tU = t U + ∆U.(6.149)6. Если k = kmax , то останавливаем процесс. Если k < kmax , то присваиваем k := k + 1 ипереходим к п. 1 второго этапа.6.2.3.Решение задач динамического деформирования упругопластических тел пошаговым интегрированиемСначала рассмотрим неявную схему интегрирования уравнений движения (6.115) сиспользованием метода Ньюмарка.
Подставляя вектор t+∆t Ü из первого равенства (6.92)в левую часть уравнения движения (6.115), получимtK̂∆U = t+∆t R̂,(6.150)134Глава 6. Основы численных методов решения задач деформирования тел . . .где введены эффективная матрица касательной жесткости K̂ и эффективный векторнесбалансированной нагрузки t+∆t R̂:tK̂ ≡ t K + a0 M,t+∆tR̂ ≡ t+∆t R + M(a1 t U̇ + a2 t Ü) − t F.(6.151)Получим уравнения движения для i-ой итерации, переписывая уравнения (6.115) сзаменой ∆U → ∆U(i) :M t+∆t Ü(i) + t+∆t K(i−1) ∆U(i) = t+∆t R − t+∆t F(i) .(6.152)Преобразуем первый член уравнения (6.152).
Из первого равенства в (6.92) получимt+∆tÜ(i) = a0 (t+∆t U(i) − t U) − a1 t U̇ − a2 t Ü.(6.153)Из (6.123) имеемt+∆tU(i) = t+∆t U(i−1) + ∆U(i) .(6.154)Подставляя выражение для вектора t+∆t U(i) из (6.154) в правую часть (6.153) и используяобозначение (6.135), получим:t+∆tÜ(i) = a0 ∆U(i) + a0 U(i−1) − a1 t U̇ − a2 t Ü.(6.155)Из первого равенства в (6.92) имеемt+∆tÜ(i−1) = a0 U(i−1) − a1 t U̇ − a2 t Ü.(6.156)Из (6.155), (6.156) получимt+∆tÜ(i) = a0 ∆U(i) + t+∆t Ü(i−1) .(6.157)Подставляя выражение для вектора t+∆t Ü(i) из (6.157) в левую часть (6.152), получимсистему линейных алгебраических уравнений для определения вектора ∆U(i) :t+∆tK̂(i−1) ∆U(i) = t+∆t R̂(i−1) ,где введены эффективная матрица касательной жесткостивектор несбалансированной нагрузки R̂(i−1) на итерациях:t+∆tK̂(i−1) ≡ t+∆t K(i−1) + a0 M,t+∆t(6.158)t+∆tK̂(i−1) и эффективныйR̂(i−1) ≡ t+∆t R − M t+∆t Ü(i−1) − t+∆t F(i−1) .(6.159)Приведем алгоритм решения динамических задач теории пластичности по неявнойсхеме интегрирования уравнений движения с использованием метода Ньюмарка.
Этап начальных вычислений состоит из трех шагов:1. Определение матрицы масс M.2. Присвоение начальных значений векторам 0 U, 0 U̇, 0 Ü.3. Выбор шага по времени ∆t, параметров α, δ метода Ньюмарка, вычисление на их основеконстант a0 , a1 , a2 , a3 , a4 .Этап пошагового интегрирования в цикле (начальное значение параметра цикла k = 1)состоит из шести шагов:1. Определение эффективной матрицы касательной жесткости t K̂ по формуле в (6.151) иразложение ее на множители t K̂ = LDLT .2. Определение эффективного вектора несбалансированной нагрузки t+∆t R̂ по формуле в(6.151).3. Определение вектора приращений перемещений узловых точек ансамбля из уравнения(LDLT )∆U = t+∆t R̂.6.2. Применение метода конечных элементов к решению задач теории .
. .1354. Проведение итерационной процедуры уточнения решения так же, как в квазистатическом анализе, с заменами матрицы касательной жесткости t+∆t K(i−1) на эффективную матрицу касательной жесткости t+∆t K̂(i−1) и вектора несбалансированной нагрузкиt+∆tR − t+∆t F(i−1) на эффективный вектор несбалансированной нагрузки t+∆t R̂(i−1) .5. Определение новых векторов ускорений, скоростей и перемещений узловых точек ансамбля по формулам (6.92) с учетом равенства t+∆t U = t U + ∆U.6. Если k = kmax , останавливаем процесс.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
