Применение численных методов для анализа процессов и объектов ОМД

Задачи механики сплошных сред сводятся к дифференциальным уравнениям в частных производных, которые необходимо интегрировать при определенных краевых условиях [12, 13, 16, 17]. Приближенное решение краевых задач во многих случаях удается получить с применением так называемых прямых методов [29]. Прямыми называются методы приближенного решения задач теории дифференциальных и интегральных уравнений, которые сводят эти задачи к конечным системам алгебраических уравнений. В теории и практике применения прямых методов особое место занимают два метода: метод Ритца и метод Галеркина.

В первом из них задача интегрирования дифференциального уравнения заменяется некоторой равносильной вариационной задачей. Второй основан на ортогонализации невязки операторного уравнения по отношению к координатной системе функций и, вообще говоря, не связан с какой либо вариационной задачей.

Достаточно подробно данные методы и применение их для решения задач ОМД рассмотрены в работах [13, 15, 16, 17, 29].

Проекционные методы решения задач ОМД

Метод Ритца

Пусть требуется найти минимум некоторого функционала J(x) с областью определения D_J.

Выберем координатную систему функций φ₁, φ₂, …, φ_n, удовлетворяющую следующим требованиям [29]:

1) элементы координатной системы, взятые в любом конечном количестве, линейно независимы;

2) координатная система полна в некоторой метрике, определенной на области D_J;

3) при любых значениях постоянных а₁, а₂,…, а_n элемент

 

(4.1)

 

принадлежит D_J и выражение J(x_n) имеет смысл.

Рассматривая его как функцию конечного числа переменных а₁, а₂,…, а_n, найдем те значения, при которых J(x_n) достигает минимума. С этой целью необходимо решить следующую систему уравнений

 

(4.2)

(необходимое условие экстремума J(x_n)). Убедившись, что найденные значения постоянных а_i действительно реализуют минимум величины J, подставим эти значения в выражение (4.1). В результате получим элемент x_n, который назовем n –м приближением по Ритцу решения данной вариационной задачи.

Для неоднородных граничных условий можно искать n –е приближение по Ритцу в следующем виде

 

 

(4.3)

где элемент φ₀ удовлетворяет заданным граничным условиям, а φ_i удовлетворяет соответствующим однородным граничным условиям.

Решение системы уравнений (4.2) является в общем случае весьма сложной задачей. Она существенно упрощается, если J(x_n) – квадратичный функционал, в этом случае уравнения (4.2) линейны относительно а_i.

На практике во многих случаях приходится ограничиваться сравнительно небольшим числом членов рядов (4.1) и (4.2), поэтому удачный выбор координатных функций имеет решающее значение. При решении вариационных задач обработки металлов давлением для выбора координатных функций обычно используют результаты экспериментальных исследований.

Пример. Рассмотрим расчет деформированного состояния полосы прямоугольного сечения при кузнечной протяжке (рис. 2.9) при указанных там граничных условиях. Модель построим для поперечного сечения yOz.

Кривую упрочнения Т(Н) аппроксимируем следующей функцией [16]

 

Т = 1,88 Н ^1/3.

 

Эта зависимость соответствует деформации стали марки 45 при 1100⁰С.

Функционал для рассматриваемого случая

 

(4.4)

 

где

 

 

Подходящей последовательностью функций вида (4.1) для поля скоростей, удовлетворяющей граничным условиям, будет

 

(4.5)

 

Ограничим (4.5) двумя членами ряда. Второй компонент скорости найдем из условия несжимаемости при

Очевидно, что из выражения (4.5) следует

 

 

После интегрирования данного выражения для нахождения и преобразований получили компоненты тензора скорости деформации [16]

 

(4.6а)

(4.6б)

 

Деформированное состояние описывается приближенно формулами (4.6), которые содержат один варьируемый коэффициент а₁, который определим из условия экстремума функционала (4.4)

 

 

В результате найдено а₁ =0,73 v, подставив которое в (4.6) получим

 

 

Распределение интенсивности скоростей деформации в безразмерном виде Н¹ = Н/(2v/h) для одной четверти представлено на рис. 4.1 (подготовлено с применением пакета EXCEL).

 

 

Шкала H1

z/h

y/b

H1

 

Рис. 4.1. Распределение интенсивностей скорости деформации

по сечению заготовки при протяжке

 

Метод Галеркина

Пусть требуется найти решение уравнения

 

 

Как и в методе Ритца, выбираем координатную систему функций φ₁, φ₂, …, φ_n, удовлетворяющую следующим требованиям [29]:

4) элементы координатной системы, взятые в любом конечном количестве, линейно независимы;

5) координатная система полна в некоторой метрике, определенной на области D_L;

6) при любых значениях постоянных а₁, а₂,…, а_n элемент

 

(4.7)

 

принадлежит D_L и выражение Lx_n имеет смысл.

Запишем условие ортогональности невязки уравнения к первым n координатным функциям

 

 

Это означает, что выражение «равно нулю» в подпространстве Н⁽ⁿ⁾ с базисом φ_i, i=1,2,…,n, т.е. ортогонально базисным функциям и любому элементу этого подпространства.

В результате получаем систему из n уравненийдля нахождения коэффициентов а₁, а₂,…, а_n. Если оператор L – линеен, то эта система представляет собой систему линейных алгебраических уравнений относительно а_i.

Решив полученную систему и подставив коэффициенты а_i в (4.7), получаем элемент x_n, который назовем n – м приближением по Галеркину решения данной задачи.

Методы Ритца и Галеркина широко используются в методе конечных элементов.

 

Метод конечных элементов

 

Среди численных методов решения наибольшее распространение получили метод конечных разностей (МКР) и метод конечных элементов (МКЭ). Поскольку для решения связанных краевых задач ОМД в настоящее время широко используется МКЭ и имеются специализированные пакеты программ, то краткое рассмотрение и решение тепловых задач МКР будет рассмотрено на практических занятиях.

 

4.2.1. Понятие о методах конечных элементов

 

Метод конечных элементов (МКЭ) в настоящее время, пожалуй, самый распространенный в мире численный метод. К его достоинствам относятся:

- возможность счета на неравномерных сетках, в двумерном и трехмерном случаях для областей сложной геометрии;

- "технологичность" методов.

Основная идея метода конечных элементов, базирующая на методах Бубнова, Галеркина и Ритца, была предложена Р.Курантом в 1943 г., но осталась незамеченной, опередив потребности практики. В 50 - х годах прошлого века с появлением первых компьютеров возникла необходимость в разработке новых инженерных подходов к численному решению задач со сложной геометрией, в которых области интегрирования разбивались на подобласти. Такие подобласти (носители финитныхбазисных функций, об этом ниже) и получили название конечных элементов.

Метод конечных элементов основан на идее аппроксимации непрерывной функции φ^(е) (в физической интерпретации - температуры, давления, перемещения, скорости и т.д.) дискретной моделью, которая строится на множестве кусочно-непрерывных функций, определенных на конечном числе подобластей, называемых конечными элементами [29, 30-32]. Исследуемая геометрическая область разбивается на элементы таким образом, чтобы на каждом из них неизвестная функция аппроксимировалась пробной функцией (как правило, полиномом). Причем эти пробные функции должны удовлетворять граничным условиям непрерывности, совпадающим с граничными условиями, налагаемыми самой задачей. Выбор для каждого элемента аппроксимирующей функции будет определять соответствующий тип элемента.

Будем рассматривать вычислительный алгоритм метода конечных элементов в формулировке, основанной на процедуре минимизации функционала, соответствующего решаемой непрерывной задаче. В результате выполнения указанной процедуры происходит замещение уравнения или системы уравнений в частных производных системой недифференциальных уравнений, имеющих в качестве коэффициентов аппроксимирующие функции, которые фактически являются значениями искомой функции в вершинах разбиения.

При решении задачи с помощью МКЭ необходимо определиться с формой конечного элемента.

Форма конечного элемента – его внешний вид, определяет точность аппроксимации границ исследуемого объекта

В одномерном случае выбор ограничен отрезком прямой. В двумерном случае форма конечного элемента может быть любой, при условии, что с помощью этого конечного элемента можно, с некоторой степенью точности аппроксимации границ, покрыть площадь произвольной формы (без перекрытия элементов). Наиболее простыми элементами для плоского случая являются треугольный и прямоугольный (со строронами, параллельными осям координат) элементы.

Для трехмерного случая форма элемента должны быть такой, чтобы с его помощью можно было бы покрыть объем произвольной формы, аппроксимировав при этом границы объекта. Наиболее простыми элементами являются тетраэдр и параллелепипед со стронами, параллельными осям координат. Условие параллельности упрощает вычисление локальных матрицы жесткости вектора нагрузок.

В общем случае алгоритм МКЭ состоит из четырех этапов:

Этап 1. Выделение конечных элементов (разбиение области на конечные элементы).

Этап 2. Определение аппроксимирующей функции для каждого элемента (определение функции элемента). На данном этапе значение непрерывной функции φ^(е) в произвольной точке е -го конечного элемента аппроксимируется полиномом

 

φ^(е)=А^(е)R + A₀, (4.8)

 

где А^(е) – вектор-строка коэффициентов полинома, A₀ – свободный член, R=(х,у,z) – вектор координат в рассматриваемой точке.

Задача этапа далее заключается в определении неизвестного вектора А^(е) и свободного члена A₀. Для этого, использя условие непрерывности функции в узлах, коэффициенты полинома выражают через вектор Ф ^(е) узловых значений функции и координаты узлов и, проделав эквивалентные преобразования, получают

 

φ^(е)= N^(e) Ф ^(е), (4.9)

 

где N^(e) – матрица–сторка, элементы которой называют функциями формы конечного элемента.

Функции Фомы легко вычисляются через координаты самой точки и координаты узлов элементов.

Этап 3. Объединение конечных элементов в анасамбль. На этом этапе уравнения (4.9), относящиеся к отдельным элементам, объединяются в ансамбль, т.е. в систему алгебраических уравнений

 

φ = NФ. (4.10)

 

Система (4.10) является моделью искомой непрерывной функции.

Этап 4. Определение вектора узловых значений функции. В общем случае вектор Ф в (4.10) вначале неизвестен. Его определние – наиболее сложная процедура в МКЭ.

Оазработано несколько алгоритмов вычисления вектора Ф. Один из алгоритмов основан на минимизации функционала, связанного с физическим смыслом решаемой задачи, и состоит из следующих этапов:

Этап 1. Выбор функционала, зависящего для стационарных задач от искомой функции φ и её частных производных по вектору пространственных координат

 

 

где V – объем.

Функционал J представляется суммой соответствующих функционалов, относящихся к отдельным конечным элементам

 

(4.11)

 

где N – число элементов.

Этап 2. Подстановка аппроксимирующего выражения (4.9) в (4.11) и вычисление производных по формулам вида

 

 

Этап 3. Минимизация по вектору Ф функционала J. Для этого составляяются уравнения

 

(4.12)

 

Суммирование выражений (4.12) по конечным элементам приводит к системе алгебраических уравнений

 

КФ=F, (4.13)

 

где К – матрица коэффициентов – матрица жесткости, F – вектор нагрузки.

Матрица жесткости и вектор нагрузки представлят математическую модель в МКЭ.

Этап 4. Решение системы (4.13), позволяющее определить неизвестный вектор узловых значений.

Найденные значения вектора Ф подставляют в (4.10), после чего значение функции φ легко вычисляются в любой точке заданной области. Процедуру определения аппроксимирующей функции элементов выполняют один раз для типичного элемента области безотносительно к его топологическому положению в ней.

Представим, что заготовка в состоянии плоской деформации разделена на конечное число треугольных призматических элементов с воображаемыми границами (рис. 4.2). Элементы соединены только в узловых точках (узлах) и силы не могут передаваться через боковую поверхность элементов. Элементы и узлы так пронумерованы, что соседние элементы или узлы имеют близкие номера.

При небольшом перемещении инструмента узлы пермещаются в новые положения. Перемещения узлов считаются неизвестными параметрами, которые предстоит опрелеоить при заданных граничных условиях. В двумерной модели неизвестных параметров вдвое больше, чем узлов, поскольку каждому узлу соответствует две компоненты перемещения u вдоль оси x и v вдоль оси y. Приращения напряжений и деформаций могут быть определены в результате вычисления перемещений.

В качестве примера рассмотрим треугольный элнмент ijm в условиях плоской деформации (рис. 4.3).

 

 

Рис. 4.2. Разбивка на конечные элементы и индексация (нумерация) узлов и элементов

 

Рис. 4.3. Треугольный элемент с тремя узлами

 

Компоненты перемещения или скорости элемента интерполируются функцией перемещения. Простейшими функциями перемещения являются линейные полиномы

 

(4.14а)

(4.14б)

 

Коэффициенты α₁, α₂, …, α₆ определяются путем подстановки в эти уравнения перемещений и координат узловых точек, т.е.

 

(4.15)

 

Решение этих уравнений позволяет найти

 

(4.16)

 

где А – площадь треугольника ijm, т.е

 

 

а a_i, b_i и c_i определены следующим выражением

 

 

Другие коэффициенты получаются путем циклической перестановки индексов. Используя эти значения, уравнения (4.14) можно переписать в терминах скоростей узлов

 

(4.17)

 

Зависимости (4.17) можно представить в следующем виде

 

(4.18)

 

или в матричном виде

 

,

 

где

 

называют функциями формы элемента.

Рассмотрим процедуру составления ансамбля конечных элементов. В соответствии с рис. 4.2 начнем с узла 1 и ведем отсчет против часовой стрелки. Соответствие между этими обозначениями и глобальными номерами узлов следующее:

элемент 1 i = 1; j = 5; m =6;

элемент 2 i = 1; j = 6; m =2;

элемент 3 i = 2; j = 6; m =7;

элемент 4 i = 2; j = 7; m =3 и т.д.

Подставляя полученные значения в (4.18), получим

 

 

Данная система является моделью искомой непрерывной функции.

Матрица жесткости элемента в (4.13) находится из условия баланса работы внешних сил и суммарной внутренней работы [30-32] и имеет вид

 

,

где V – объем элемента, а матрица В определяется через производные от функций (4.17)и имеет вид

 

 

Матрица D вытекает из соотношения между напряжениями и деформациями σ = D ε и при плоском напряженном состоянии для условия пластичности Мизеса имеет вид

 

 

где Е – модуль упругости Юнга, v – коэффициент Пуассона, G = E/2(1+v) –модуль сдвига, – компоненты девиатора напряжений, – коэффициент упрочнения , – напряжение текучести.

Общая матрица жесткости для ансамбля элементов выразится в виде

 

 

и будет представлять собой симметричную (в силу теоремы Бетти) матрицу. В итоге решение для (4.13) сведётся к решению системы уравнений с разряженной (ленточной)матрицей.

Мы рассмотрели только основы теории метода конечных элементов, для более глубокого изучения особенностей применения метода для решения задач ОМД необходимо обратиться к работам [29-32], а также уточнить эти вопросы при практическом освоении пакта ANSYS.

Пример решения одномерной задачи с помощью МКЭ.

Пусть необходимо найти удлинение балки, с одним закрепленным концом (рис. 4.4) с продольной нагружающей силой.

 

Рис. 4.4. Схема балки с одним закрепленным концом и продольной нагружающей силой

Уравнение, описывающее состояние балки имеет вид:

здесь y - удлинение, F - нагружающая сила, S - площадь поперечного сечения, E - модуль Юнга.

В соответствии с алгоритмом решения задач с помощью МКЭ:

1. Выбираем конечный элемент. Для одномерной задачи выбор ограничен только отрезком прямой.

2. Выбираем функцию формы конечного элемента, то есть фактически выбираем аппроксимацию решения внутри конечного элемента. Будем считать, что удлинение внутри конечного элемента меняется по линейному закону:

 

(4.19)

 

Предполагаем, что нам известны узловые значения удлинений Y_i и Y_j (см. рис. 4.5).

Из (4.19) при , а при .

Из данной системы уравнений находим значения и и подставляем в (4.19), выделяя коэффициенты при и :

 

Рис. 4.5. Схема узловых значений удлинений

 

,

где - вектор функции формы конечного элемента, его составляющие элементы - глобальные базисные функции, отличные от нуля в пределах этого элемента.

3. Разбиваем область на конечные элементы.

4. Получение локальных матрицы жесткости и вектора нагрузок конечного элемента.

Локальная матрица жесткости и вектор нагрузок - математическая модель конечного элемента. Эти термины употребляются не только в задачах строительной механики, но и в других предметных областях

Фактически для их получения необходимо применить метод взвешенных невязок в пределах конечного элемента с аппроксимацией, полученной в в соответствии с методом Галеркина:

 

Раскрываем интеграл в предположении, что площадь поперечного сечения элемента постоянна:

 

 

Приводим уравнение к следующему виду:

(4.20)

 

Получили локальные матрицу жесткости и вектор нагрузок.

5. Ансамблирование.

Интеграл по одному конечному элементу мы вычислили в (4.20).

Глобальная матрица жесткости будет иметь размерность, определяемую числом узлов сетки, в нашем примере – 4. Вектор неизвестных составляют перемещения в этих узлах. Локальная матрица жесткости каждого конечного элемента даст аддитивный вклад в глобальную матрицу в соответствии с узлами подключения конечного элемента (это же касается и вектора нагрузок)

 

6. Учет граничных условий. В нашем примере , то есть можно вычеркнуть первый столбец и первую строку.

7. Решение системы уравнений

 

 

В результате найдем удлинение в каждом узле.

 

 

⇐ Предыдущая45678910111213Следующая ⇒

Поделиться:

Воспользуйтесь поиском по сайту: