Учебное пособие по курсу лекций (1164067), страница 58
Текст из файла (страница 58)
ДлявыводаопределяющихсоотношенийМКЭвариационнымметодом377используют матричную запись уравнений. Математическая модельконечного элемента для задач механики деформируемого тела в этом случаедолжна представлять зависимость напряжений {σ } в элементе от значенийузловых неизвестных. Такую модель следует получить в матричном виде.{σ } =f({ u }) или {σ } = f ({ v })eeВ общем виде такую задачу можно решить, используя уравнения Кошии физические уравнения связи напряженного и деформированногосостояний. Мы до настоящего момента использовали тензорную запись такихуравнений.Уравнения Коши для определения деформаций (скоростейдеформаций) по полю перемещений (скоростей) материальных частиц втензорном виде:11ε ij = ui, j + u j ,i ,ε ij = vi, j + v j ,i ,i , j = x, y , z22В матричной записи уравнения Коши для деформаций выглядятследующим образом:{ε } = [∂ ]{u} ,где {ε } - вектор-столбец, составленный из компонентов тензорадеформаций, [∂ ] - матрица дифференциального оператора.
Например, дляплоской задачи εx=∂ux/∂x, εy=∂uy/∂y, γxy=∂ux/∂y+∂uy/∂x, поэтому:⎡∂⎤0⎥⎢⎧ ε x ⎫ ⎢ ∂x⎥∂ ⎥ ⎧u x ⎫⎪ ⎪ ⎢×⎨ ⎬⎨εy ⎬= ⎢ 0⎥∂y⎩u y ⎭⎪γ ⎪ ⎢⎩ xy ⎭∂ ∂ ⎥⎢⎥⎣ ∂ y ∂x ⎦Для осесимметричной задачи:⎡∂⎤⎢ ∂ρ 0 ⎥⎥⎧ερ ⎫ ⎢∂ ⎥⎢⎪ε ⎪0⎪ z ⎪ ⎢∂z ⎥ × ⎧u ρ ⎫=⎨ε ⎬ ⎢ 1⎥ ⎨ ⎬0 ⎥ ⎩u z ⎭⎪ θ ⎪ ⎢⎪⎩γ ρz ⎪⎭ ⎢ ρ⎥∂ ⎥⎢∂⎢⎣ ∂z ∂ρ ⎥⎦Используя функции формы для аппроксимации поля перемещенийконечного элемента получим:()(){ε } = [∂ ]{u} = [∂ ][ N ]{u (e) } = [ B ]{u (e)} ,378где [B] - матрица градиентов - частных производных функций формы попространственным координатам.
Например, для 2D треугольного линейного элементаматрица градиентов:⎡∂⎢⎢∂ x⎢[ B] = ⎢ 0⎢⎢∂⎢∂ y⎣⎤0 ⎥⎥∂ ⎥ ⎡ Ni×⎢∂ y ⎥⎥ ⎢⎣ 0∂ ⎥∂ x ⎥⎦0Nj0NkNi0Nj00 ⎤⎥=N k ⎥⎦⎡bi 0 b j 0 bk 0 ⎤⎥1 ⎢=⎢ 0 ci 0 c j 0 ck ⎥2∆ ⎢⎥⎣⎢ci bi c j b j ck bk ⎦⎥Уравнения Коши для скоростей деформаций имеют следующий вид:{ε } = [∂ ]{v} = [∂ ][ N ]{v(e)} = [ B ]{v(e)}Здесь матрица [B] имеет тот же смысл.Уравнения связи напряженного и деформированного состояний(физические уравнения) для упругих деформаций (обобщенный закон Гука) втензорной записи:sijeij =, i , j = x, y , z2Geij = ε ij + δ ij ε cp ,sij = σ ij + δ ijσ cpε cp =σ cp,1ε cp = ε ii ,31σ cp = σ ii33K(Здесь sij - девиатор напряжений, eij - девиатор деформаций.
)Эти уравнения можно разрешить относительно компонентов тензоранапряжений и представить в матричном виде:{σ } = [D]{ε }Здесь {σ } - вектор столбец, составленный из компонентов тензоранапряжений, [D ] - матрица упругих коэффициентов. В частности дляплоского деформированного состояния в развернутом виде эту запись можнопредставить так77:77При плоском деформированном состоянии ε z = 0;(σz = µ σx +σ y)379µ⎡ 1⎤01− µ⎢⎥⎧ ε x ⎫⎧σ x ⎫E (1 − µ ) ⎢ µ⎥⎪ ⎪⎪ ⎪10⎨σ y ⎬ =⎢ 1− µ⎥⎨ ε y ⎬()()112µµ+−⎪τ ⎪⎢⎥ ⎪γ xy ⎪()12µ−⎩ ⎭⎩ xy ⎭0⎢ 02(1 − µ )⎥⎦⎣Преобразуя матричную запись обобщенного закона Гука с учетомранее полученных соотношений для деформаций, окончательно получимследующую матричную форму зависимости напряжений от узловыхперемещений:{σ } = [D][B] u (e)Для пластических деформаций (при пренебрежении упругими)физические уравнения приобретают вид уравнений течения (Леви-Мизеса).
Втензорном виде :3 εiε ij =sij2σi(Здесь ε ij - тензор скоростей деформаций, ε i - интенсивность скоростей{ }деформаций, σ i - интенсивность напряжений.)Матричная запись уравнений течения имеет вид:2σ{s} = i ⎡⎣ D p ⎤⎦ {ε }3ε iМатрица связи [Dp] постоянна и зависит от вида напряженногосостояния. Например, для плоского деформированного состояния.⎡1 0 0 ⎤⎡ D p ⎤ = ⎢0 1 0 ⎥⎣ ⎦ ⎢⎥⎢⎣0 0 0.5⎥⎦Окончательно, с учетом матричной записи скоростей деформаций{s} = ⎡⎣ D p ⎤⎦ [ B ] v(e)Таким образом, получено матричное уравнение, связывающеекомпоненты тензора напряжений (для упругих деформаций) или девиаторанапряжений (для вязко-пластической задачи) в любой точке конечногоэлемента в зависимости от узловых неизвестных (соответственноперемещений или скоростей).
Для треугольного симплекс элемента матрицаупругих коэффициентов [D] и матрица градиентов [B] – постоянны, поэтомунапряжения в треугольном симплекс элементе для упругих деформацийпостоянны. Это надо иметь в виду при анализе полученных результатов.{ }4.9.8 Вариационный метод построения разрешающих уравнений.Упругие деформации.Ранее был использован прямой метод построения разрешающейсистемы уравнений МКЭ для упругого деформирования. Получим этусистему с использованием вариационного метода. В вариационном методе380задачу ищут путем минимизации некоторого функционала78, связанного сфизической сущностью задачи. При построении разрешающих уравненийМКЭ для задач теории упругости используют функционал полной энергии.Π = Λ − Α,где П – полная энергия, Λ - энергия деформации (работа внутренних сил), A - работавнешних сил.Согласно вариационному принципу Лагранжа (принципу минимумапотенциальной энергии) полная потенциальная энергия тела, находящегося вравновесии под действием внешних сил, достигает минимума на истинномполе перемещений.∂Π=0∂uПри дискретизации МКЭ истинное поле перемещений замененосовокупностью перемещений узлов конечных элементов.
Поэтому вдискретном аналоге вариационного принципа Лагранжа дифференцированиевыполняют по узловым перемещениям∂Π=0∂ ueПолная энергия в дискретном случае вычисляется суммированиемсоставляющих по конечным элементам:{ }Π = ∑ Λ e − ∑ Αeeeeгде Λ - энергия деформации одного КЭ, Ae - работа внешних сил дляодного КЭРабота упругой деформации в тензорном виде:1Λ = ∫ σ ij ε ij dV2VВ матричной форме1TΛ = ∫ {ε } {σ }dV2VДействительно, например для двумерного случая:⎧σ x ⎫⎪ ⎪T∫ σ ij ε ij dV = ∫ (σ xε x + σ yε y + τ xyγ xy )dV = ∫ {ε x , ε y ,γ xy }⎨σ y ⎬dV = ∫ {ε } {σ }dVVVVV⎪σ ⎪⎩ z⎭Таким образом, работа деформации одного элемента (интегрированиеведется по объему элемента Ve):78Функционалом Ф, зависящем от функции f называется такая переменнаявеличина, которая принимает конкретное числовой значение, приподстановке в нее каждой функции f из некоторого класса функций.3811{ε }Τ {σ }dV(e) 2Λe = ∫VРанее были получены матричные выражения, связывающиедеформацию и напряжение в произвольной точке элемента с узловымиперемещениями:{ε } = [B ] u (e) , {σ } = [D ][B ] u (e)Используя правила транспонирования матричных произведений,получим79:{ }{ε }T{ }= u (e)T{ }[ B ]TТогда выражение для энергии деформации элементаследующий вид80:T11TTΛe = ∫ {ε } {σ }dV = ∫ u e [ B ] [ D ][ B ] u e dV =2 e2 eV={ }1 e Tu2V{ }T∫ [ B ] [ D ][ B ] dVVeпримет{ }{ } { }ue ={ }1 e T ⎡ e⎤ euk u⎣ ⎦2⎡k e ⎤⎣ ⎦Выражение в угловых скобках справедливо для вычисления матрицжесткости конечных элементов любой размерности и степенейаппроксимации поля перемещений в упругих задачах.
Таким образом,матрица жесткости конечного элемента в общем случае производитсяинтегрированием матричного соотношения по объему элемента:⎡ k e ⎤ = ∫ [ B ]T [ D ][ B ] dV⎣ ⎦eVВ частном случае, когда матрицы градиентов и упругихкоэффициентов являются постоянными, их можно вынести за знак интегралаи выражение упрощается:⎡ k e ⎤ = [ B ]T [ D ][ B ]V e⎣ ⎦79Здесь использовано:TЗнак [ ] после матрицы означает ее транспонирование, т.е. заменастолбцов матрицы ее строками. После транспонирования матрица [ A]размером m × n превращается в матрицу [ A] размером n × m , причем:TaijT = a jiДля транспонированных матриц справедливо следующее соотношение:([A][B ])T = [B]T [A]T80Вектор узловых перемещений вынесен из-под знака интеграла, посколькуон не зависит от координат.382Для более сложных КЭ выполняют численное интегрирование.Следует обратить внимание, что в выражении для матрицы жесткостиконечного элемента присутствует его объем.
Именно поэтому при заданиисвойств одномерных и плоских конечных элементов приходится вводитьдополнительные параметры (площадь и толщину элемента соответственно),позволяющие вычислить объем элемента.Работа внешних сил в тензорной записи имеет вид:A = ∫ piui dF ,Fpгде Fp – часть внешней поверхность на которой заданы внешние силы, pi – удельныесилы на внешней поверхности, в проекциях на координатные оси, ui – перемещенияматериальных точек внешней поверхности в проекциях на координатные оси.Конечные элементы взаимодействуют с внешней средой только вузлах, поэтому интегрирование по внешней поверхности можно заменитьсуммированием произведения узловых внешних сил на соответствующиеузловые перемещения.
Например, для двумерного треугольного конечногоэлемента с узлами i,j,k:Αe = Pxiu xi + Pyiu yi + Pxj u xj + Pyj u yj + Pxk u xk + Pyk u yk =⎧ Pxi ⎫⎪P ⎪⎪ yi ⎪⎪P ⎪⎪ xj ⎪(e) T= ⎡⎣u xi , u yi , u xj , u yj ,u xk , u yk ⎤⎦ ⎨=uP ( e)⎬⎪ Pyj ⎪⎪P ⎪⎪ xk ⎪⎪⎩ Pyk ⎪⎭Подставив их в уравнение для функционала полной энергии, получимвыражение для полной энергии системы, составленной из совокупностиконечных элементов:{ }{ }383{ } [ ]{ }{ }{ }T1 eT e eu k u − ∑ u e Pee 2eПроанализируем физический смысл полученного выражения.
Реальнаядеформирующаяся упруго система была заменена совокупностью конечныхэлементов, взаимодействующих друг с другом в узлах. Свойства каждогоэлемента определены его матрицей жесткости. Таким образом, исходнаясистема с бесконечным числом степеней свободы заменена дискретнойсистемой, неизвестными в которой являются узловые перемещения. Какойфизический принцип положить в основу определения этих неизвестныхперемещений? Естественно потребовать, чтобы сконструированное намидискретное поле перемещений максимально соответствовало реальномуполю перемещений.Для этих целей удобно использовать дискретный вариантвариационного принципа Лагранжа, рассмотренный нами ранее. Согласноэтому принципу наилучшее приближение к истинному полю перемещенийдостигается, если частная производная от полной энергии по каждомуузловому перемещению становится равной нулю:T1 e T ⎡ e⎤ e∂Π∂uk u − ∑ uePe ==∑⎣ ⎦∂ ue∂ ue e 2eΠ=∑{ } { }={ }{ } { }{ }{ }{ }{ }{ }1 e T ⎡ e⎤ e∂e Tuku−uPe = 0∑∑ee⎣⎦∂ u e 2∂ u e∂{ }{ }Воспользовавшисьправиламидифференцированияматричных81произведений , получим:⎞∂Π ⎛= ⎜⎜ ∑ k e u e ⎟⎟ − ∑ P e = [K ]{U } − {R} = 0e∂u⎝ e⎠ eИли окончательно определяющие соотношения МКЭ для упругойзадачи:{R} = [ K ]{U } ,что тождественно полученным ранее прямым методом.Здесь {R} - вектор столбец узловых нагрузок размерностью N×1,представляет собой сумму векторов внешних сил, приложенных к узламсетки.
N – число узловых неизвестных. Обычно N=r×M, где r – размерностьзадачи (r=1 для 1D, r=2 для 2D, r=3 для 3D), M – суммарное количество узлов.[K] - глобальнаяматрицажесткостисистемы,полученнаяпоэлементным суммированием матриц жесткости каждого элемента82[ ]{ }{ }(){ }()∂∂A TBA = 2BA ,A TB = B∂A∂A82Поэлементное суммирование, или ансамблирование выполняют поспециальному алгоритму, который будет рассмотрен иже.81384[ K ] = ∑ [k ( e ) ] =∑ ∫ [ B]Τ [ D][ B]dVee VeПоскольку глобальная матрица жесткости получается суммированиемэлементов матриц жесткостей элементов, то она наследует их свойства(положительную определенность и симметричность) и кроме того имеетленточную структуру - т.е.