1625914689-e957c8b7a8e4003fe3539e4e0e465a65 (532400), страница 31
Текст из файла (страница 31)
. . a4 по формулам (6.91).4. Определение эффективной матрицы жесткости K̂ по выражению в (6.94).5. Выполнение триангуляризации эффективной матрицы жесткости K̂ = LDLT .На втором этапе проводятся вычисления в циклах, соответствующих шагам повремени. На каждом шаге выполняются три операции:1. Определение эффективного вектора внешних сил2. Определение вектора перемещенийt+∆tt+∆tU из системы алгебраических уравнений:(LDLT ) t+∆t U = t+∆t R̂.3. Определение векторов t+∆t U̇ ипо времени t → t + ∆t.t+∆tR̂ по выражению в (6.94).(6.95)Ü из формул (6.92), переход к следующему шагу126Глава 6.
Основы численных методов решения задач деформирования тел . . .6.2.Применение метода конечных элементов к решению задач теории пластичности6.2.1.Пространственная дискретизация уравнений квазистатического/динамического движения упругопластических телЗапишем уравнение принципа виртуальных работ в момент времени t + ∆t:∫t+∆tσij δεij dV = t+∆t R̂.(6.96)VВведем обозначения:∆σij ≡ t+∆t σij − t σij ,∆εij ≡ t+∆t εij − t εij ,∆ui ≡ t+∆t ui − t ui .(6.97)Из (6.97) получим следующие равенства:1∆εij = (∆ui,j − ∆uj,i ),2t+∆tσij = σij + t σij .(6.98)Определяющие соотношения упругопластичностиepσ̇ij = Cijklε̇kl(6.99)линеаризуем относительно момента времени t:ep∆σij = t Cijkl∆εkl .(6.100)epσij = t σij + t Cijkl∆εkl .(6.101)Из (6.98), (6.100) имеемt+∆tИз (6.96), (6.101) получим линеаризованное относительно момента времени t уравнениепринципа виртуальных работ:∫∫t ept+∆tCijkl ∆εkl δεij dV =R̂ − t σkl δεij dV.(6.102)VVВводим векторы t σ̄, ∆ε̄ и симметричную матрицу t C,ep:тензора четвертого ранга t Cijkl t tσ11∆ε11C11 t C12tt ∆ε22 σ22 C22 tσ33 tt ∆ε33 σ̄ ≡ t σ12 , ∆ε̄ ≡ 2∆ε12 , C ≡ симм. t 2∆ε13 σ13 t2∆ε23σ23составленную из компонентtC13tC23tC33tC14tC24tC34tC44tC15tC25tC35tC45tC55tC16tC26tC36tC46tC56tC66.(6.103)Перепишем уравнение (6.102) с помощью введенных в (6.103) обозначений:∫∫T tt+∆tδε̄ C∆ε̄ dV =R̂ − δε̄T t σ̄ dV.V(6.104)VЭто уравнение можно использовать для описания динамического и квазистатического деформирования тел из упругопластического материала.6.2.
Применение метода конечных элементов к решению задач теории . . .127При квазистатическом деформировании рассматриваем функционал (5.197), который в момент времени t с использованием введенных векторов и матрицы имеет следующий вид (добавим потенциал сосредоточенных сил):∫J(u̇) = ε̄˙ T t Cε̄˙ dV − Ř,(6.105)Vгде(6.106)Ř = ŘV + Řp + Řc .Здесь∫∫Tu̇ ḃρ dV,ŘV =u̇ ṗ ρ dS,Řp =V∗TŘc =Kc∑Ṙik u̇ki ,(6.107)k=1Spгде векторы u, b и p введены в (6.1).Конечно-элементную аппроксимацию функций, входящих в функционал J(u̇), проводим точно так же, как это делали при аппроксимации величин, входящих в правую частьвыражения (6.21), для функционала I(u) (см.
п. 6.1.1). Действуя по схеме, представленной в п. 6.1.1, с заменой U → U̇, получим систему обыкновенных дифференциальныхуравнений:tKU̇ = Ṙ,(6.108)где введена матрица касательной жесткости ансамбля узловых точек t K, которая находится суммированием матриц касательных жесткостей элементов и вектора скоростивнешних сил ансамбля узловых точек видаK∑˙ j + R̃˙ j ) + Ṙ .Ṙ =(R̃cpV(6.109)j=1Модифицированные матрицы касательной жесткости элементов t K̃e и векторы скоростей˙ e и R̃˙ e находятся из матриц касательных жесткостей элементоввнешних сил элементов R̃pVt eK и векторов скоростей внешних сил элементов ṘeV и Ṙep с помощью булевых матрицAe по выражениям, аналогичным выражениям, представленным в (6.31).
В свою очередь,матрицы t Ke и векторы ṘeV , Ṙep определяются следующим образом:eT tK =∫∫∫tBṘeVCB dV,TH ḃρ dV,=ṘepHT ṗ∗ dS.(6.110)SpeVeVe=К уравнениям (6.108) добавляем начальные условия:0U = U0 при t = 0,(6.111)где U0 – заданный вектор начальных перемещений (обычно принимается U0 = 0).МКЭ можно рассматривать как специальную форму метода Бубнова – Галёркина. В этом случае перепишем линеаризованное выражение принципа виртуальных работ(6.104) в следующем виде:K ∫∑j=1VjT tδε̄C∆ε̄ dV =K∑j=1t+∆t(R̂Vj−t+∆tjR̂M+t+∆tR̂pj )+t+∆tR̂c −K ∫∑j=1Vjδε̄T t σ̄ dV. (6.112)128Глава 6. Основы численных методов решения задач деформирования тел .
. .Для каждого из слагаемых в (6.112) проводим конечно-элементную аппроксимацию.Получаем следующие аппроксимированные выражения (кроме виртуальной работы сосредоточенных сил, так как выражение для t+∆t R̂c является точным):∫∫T teT t e eδε̄ Cε̄ dV = δU K U ,δε̄T t σ̄ dV = δUe T t Fe ,(6.113)Vet+∆tR̂Vet+∆tVeet+∆t e= δURV ,R̂M = δUe T Me t+∆t Üe ,R̂pe = δUe T t+∆t Rep , t+∆t R̂c = δUT t+∆t Rc .e T t+∆tЗдесь использовали обозначения матрицы касательной жесткости элемента t Ke , введенной в (6.110), вектора внутренних сил элемента t Fe , введенного в (6.82), матрицы массэлемента Me , которая введена в (6.38), и векторов объемных и поверхностных сил элемента t+∆t ReV , t+∆t Rep , введенных в (6.26), (6.28). Отметим, что вектор сосредоточенныхсил t+∆t Rc на уровне элементов не вычисляется, а сразу определяется как вектор сосредоточенных сил ансамбля узловых точек заполнением ненулевых элементов значениямисосредоточенных сил по соответствующим степеням свободы.
В (6.113) введен также вектор приращений узловых перемещений элемента:∆Ue ≡ t+∆t Ue − t Ue .(6.114)Действуя подобно тому, как действовали при выводе системы обыкновенных дифференциальных уравнений (6.43), приходим к следующей системе обыкновенных дифференциальных уравнений:M t+∆t Ü + t K∆U = t+∆t R − t F,(6.115)где ввели вектор приращений перемещений ансамбля узловых точек:∆U ≡ t+∆t U − t U.(6.116)К системе (6.115) добавляем начальные условия (6.44).Сравнивая уравнения (6.108) и (6.115), видим, что при отбрасывании динамического члена уравнение (6.115) не сводится к уравнению (6.108).
Уравнение (6.108) можноинтегрировать любой численной процедурой (например, процедурой метода Рунге – Кутты). В частном случае, используя метод Эйлера, получим приближенные выражения длявекторов скоростей узловых точек ансамбля U̇ и скоростей внешних сил Ṙ:U̇ ≈ (t+∆t U − t U)/∆t = ∆U/∆t,Ṙ ≈ (t+∆t R − t R)/∆t.(6.117)Подставляя выражения для U̇ и Ṙ из (6.117) в (6.108), получимtK∆U = t+∆t R − t R.(6.118)Из уравнений равновесия следует, что t R = t F, поэтому из (6.118) имеемtK∆U = t+∆t R − t F.(6.119)Это уравнение можно также получить, отбрасывая первый член (соответствующий инерционным силам) в (6.115).Система обыкновенных дифференциальных уравнений (6.115) служит основой применения неявных схем ее интегрирования. Мы получили эту систему, исходя из записиуравнения принципа виртуальных работ (6.5) в момент времени t + ∆t в виде (6.96), дальнейшей линеаризацией этих уравнений относительно момента времени t в виде (6.104) и6.2.
Применение метода конечных элементов к решению задач теории . . .RR3R129приближенное решениеточное решение42R1R00 1U2U34UUUРис. 6.2. Пошаговое интегрирование уравнений движения (равновесия) без применения итерационной процедуры уточнения решения.применением последующей конечно-элементной аппроксимации входящих в эти уравнения неизвестных. Для применения явных схем интегрирования уравнений динамическогодеформирования упругопластических тел записываем уравнение (6.5) в момент времениt:∫tσij δεij dV = t R̂.(6.120)VПроводя конечно-элементную аппроксимацию выражений в (6.120), приходим к системеобыкновенных дифференциальных уравнений:M t Ü = t R − t F.(6.121)К этим уравнениям добавляются начальные условия (6.44).6.2.2.Решение задач квазистатического деформирования упругопластических тел пошаговым интегрированиемПроцедуру пошагового решения уравнений квазистатики для задачи с одной степенью свободы иллюстрирует рис. 6.2.
Недостатком такой процедуры интегрирования уравнений является то, что при относительно большом шаге ∆t приближенное решение можетдостаточно сильно отклониться от точного. Во избежание этого необходимо применятьитерационные процедуры уточнения решения.Среди методов итерационного уточнения решения наиболее известным являетсястандартный метод Ньютона – Рафсона. В этом методе на каждой итерации решаетсясистема линейных алгебраических уравнений:t+∆tK(i−1) ∆U(i) = t+∆t R − t+∆t F(i−1) .(6.122)Здесь индекс в скобках означает номер итерации, правая часть уравнения называетсявектором невязки, или вектором несбалансированной нагрузки, вектор разности перемещений узловых точек ансамбля для двух соседних итераций обозначается следующимобразом:(6.123)∆U(i) ≡ t+∆t U(i) − t+∆t U(i−1) .Начальные условия итерационного процесса имеют следующий вид:t+∆tK(0) = t K,t+∆tU(0) = t U,t+∆tF(0) = t F.(6.124)Итерационная процедура заканчивается тогда, когда выполняются некоторые критериисходимости итерационного процесса, при этом норма вектора ∆U(i) становится малой величиной.
Схема процесса уточнения решения стандартным методом Ньютона – Рафсона130Глава 6. Основы численных методов решения задач деформирования тел . . .абRt+Dtt+Dtt+Dtt+Dt(1)R- Ft+DtR(2)R- Ft+DtRttt+Dtt+Dtt+Dt(1)(1)(2)KKttFt+DtR- FRtKt+DtR- FKF000tUt+DtU(1)t+Dt(2)Ut+DtUU0tUt+DtU(1)t+Dt(2)Ut+DtUUРис. 6.3. Иллюстрация применения стандартного a) и модифицированного б ) методов Ньютона– Рафсонапоказана на рис. 6.3, а.
Недостатком итерационной процедуры стандартного метода Ньютона – Рафсона является то, что на всех итерациях при решении системы уравнений (6.122)надо формировать матрицу жесткости и проводить ее факторизацию. В модифицированном методе Ньютона – Рафсона матрица касательной жесткости на каждой итерации неперевычисляется. Вместо этого на каждой итерации используется одна и та же матрицаtK, сформированная для какого-то одного момента времени t. В частности, это можетбыть матрица 0 K.