1 марта 2003 г.

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

Здесь приводится изложение упрощенного алгоритма решения плоской задачи механики деформируемого твердого тела методом конечных элементов в пакете Mathcad, опубликованного в моей статье в журнале Exponenta.Pro (№3, 2003 г.), а также на форуме Exponenta.ru. Дату поста ставлю соответствующую.

Рассматривается простейший и в то же время наиболее распространенный вариант с разбиением области на треугольные элементы. (По ходу дела ориентировался на алгоритм, изложенный в книге Фадеев А.Б. Метод конечных элементов в геомеханике. - М.: Недра, 1987).


1. Подготовка исходных данных.

Поскольку необходимо задать информацию о каждом из элементов и узлов расчетной области, удобнее всего использовать для подготовки данных табличный редактор Excel, тем более что в Mathcad предусмотрена возможность импорта данных из файлов формата .prn. Создается в Excel два файла с таблицами, содержащими сведения об элементах и узлах. Структура таблиц и размерности величин в них приведены на рис. 1 и 2. В таблице данных об узлах имеются два столбца специальных переменных Рх и Ру, которым присваивается признак фиксации перемещения вдоль оси 0х или 0у соответственно (принимает значение 1, если задано нулевое перемещение и 0 – если перемещение неизвестно).


Рис. 1. Структура таблицы исходных данных с информацией об элементах.


Рис. 2. Структура таблицы исходных данных с информацией об узлах и заданных узловых силах и перемещениях.

Для сохранения таблиц в нужном формате, выбираем Файл->Сохранить как…, указываем в соответствующих полях имя файла и тип – Форматированный текст (разделитель – пробел). После нажатия кнопки Сохранить, нажать в появившемся окне Да. Таким образом, получаем файлы с именами, например, EL_1.prn и KN_1.prn.

2. Загрузка данных в Mathcad. Подготовка переменных.

Для удобства нумерации элементов массивов далее в книге Mathcad индекс первых элементов в массивах устанавливается равным единице:

ORIGIN:=1

Для получения данных из файлов в Mathcad используется функция READPRN(“filename.prn”) (возможно указание полного пути к файлу, в противном случае используется текущая папка, путь к которой можно узнать с помощью функции CWD).



Предположим, ранее созданные файлы находятся в папке DATA на диске D: Присваивается их содержимое матрицам DEL и DKN:



Присвоим соответствующим переменным значения из матриц: 




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

Площадь n–го элемента удобно задать в виде функции пользователя (в векторе V перечисляются глобальные номера узлов элемента):



Общая площадь расчетной области:



В массовые силы пересчитывается вес элементов, поровну на каждый из их узлов:




Узловые силы, перемещения и их признаки размещаются в векторах последовательными парами значений: на четных позициях вертикальные компоненты, на нечетных горизонтальные: 



3. Расчет матрицы жесткости системы.

Матрица жесткости системы [Kc] получается путем объединения матриц жесткости элементов [K], которые, в свою очередь рассчитываются по следующему выражению



где Delta - площадь элемента; [B] - матрица производных функций формы (функция влияния узлов), [D] - матрица связи напряжений и деформаций:



Площадь элемента вычисляется ранее заданной функцией пользователя A(n). Матрицу [D] также удобно задать в виде функции пользователя; для условий плоской деформации она будет иметь вид




Матрица [B] связывает между собой перемещения узлов элемента с его деформацией:



где







где






(выражения для функций формы Nj, Nk получаются путем круговой подстановки индексов в порядке i, j, k),

i, j, k - номера узлов элемента, xi,j,k, yi,j,k - координаты узлов.


После несложных преобразований матрицу [B] можно представить в виде



Матрицу [B] представим в виде функции пользователя, задав предварительно вспомогательную матрицу P, определяющую порядок перестановки индексов в функциях формы:





Матрица жесткости системы вычисляется в следующем программном блоке:





(объединение матриц жесткости элементов в МЖС производится по следующему правилу: член МЖС Kci,j является суммой членов Ki,j из матриц жесткости всех элементов, примыкающих к узлу с і-й степенью свободы).  


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

Основное уравнение МКЭ представляет собой систему линейных алгебраических уравнений (СЛАУ). Если какое-то і-е узловое перемещение задано, то число неизвестных уменьшается на единицу. При этом элементы і-го столбца МЖС следует умножить на заданное перемещение и результаты умножения вычесть из заданных узловых сил:



После этого і-й столбец и і-я строка МЖС, а также і-й неизвестный член в векторе сил могут быть удалены. Для удаления строк и столбцов из МЖС используем субматрицы, заданные функциями пользователя; М11 - удаляет первую строку и столбец, Mnn - последние, МI-IV - промежуточные.







Таким образом, нужно решить систему линейных алгебраических уравнений (СЛАУ) . В данном случае возможности системы Mathcad позволяют сильно упростить задачу. Для этого предусмотрена функция lsolve(M,V) для нахождения вектора решения СЛАУ, коэффициенты которой содержатся в массиве М, а свободные члены – в векторе V.



Программный модуль слева возвращает «на места» в общем векторе заданные узловые перемещения, ранее из него удалённые. Второй блок создаёт два вектора с осевыми компонентами узловых перемещений.
 

5. Нахождение осевых деформаций и напряжений в элементах

Зная полученные узловые перемещения, можно рассчитать для каждого элемента деформации и напряжения по соотношениям, упомянутым в п.3 (сигма и эпсилон):



В каждом элементе также подсчитываются главные напряжения и угол между осью 0у и вектором максимального главного напряжения . Чтобы избежать деления на ноль в строке вычисления угла использовано условное выражение, которое в случае равенства нулю выражения в знаменателе дроби присваивает углу значение .
.
 

6. Сохранение результатов.

Расчет по вышеизложенной процедуре занимает довольно непродолжительное время (например, на ПК с процессором Pentium-IV-1300 MHz; 128 MB RAM время расчета для области из 119 элементов 95 узлов составляет ~3 сек.), тем не менее, желательно сохранить результаты для последующего анализа.

Для этого сформируем матрицы, характеризующие НДС и поле перемещений, записав в них также координаты центров элементов и узлов:




(для нахождения центров элементов использована функция mean( ), возвращающая среднее значение элементов вектора)



Для записи данных в файл в Mathcad предусмотрена функция WRITEPRN(«filename.prn»); перед её использованием можно задать предварительно количество знаков после запятой переменной PRNPRECISION и ширину столбца в файле переменной PRNCOLWIDTH:



Дальнейший анализ данных наиболее удобен в графическом виде. 


7. Пример расчета.

Выполним расчет НДС вокруг выработки круглого сечения, залегающей на глубине Н от поверхности. В целях упрощения принимаем расчетную область в виде квадрата со сторонами 6*D (D - диаметр выработки); нагрузка от веса вышележащих горных пород задается в виде распределенных напряжений на границах области: , . В силу симметрии достаточно рассмотреть четверть области (рис. 3).


Рис. 3. Расчетная схема и её конечно-элементное представление.

В данном случае при разбиении на треугольные элементы получилась сеть из 95 узлов и 119 элементов. Нумерация – произвольная. 


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


В силу симметрии граничные условия по перемещениям следующие: горизонтальные компоненты вдоль вертикальной (x=0) и вертикальные вдоль горизонтальной (y=0) сторон квадрата равны нулю. Неизвестны перемещения всех узловых точек внутри массива, на контуре выработки и на грани области. 


Результаты расчета можно представлять в виде эпюр (рис. 4), изолиний напряжений или перемещений (рис. 5, а), поверхностей уровня (рис. 5, б). Сохранение и представление результатов расчета в виде векторов (матриц) позволяет сделать это без затруднений.


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






Рис. 5. Примеры визуализации результатов расчета.


Литература:
  1. Фадеев А.Б. Метод конечных элементов в геомеханике. – М.: Недра, 1987. – 221с.
  2. Ержанов Ж.С., Каримбаев Т.Д. Метод конечных элементов в задачах механики горных пород. – Алма-Ата: Наука, 1975. – 239 с.
  3. Зинкевич О. Метод конечных элементов в технике: Пер. с англ. – М.: Мир, 1975. - 542 с.
  4. Норри Д., де Фриз Ж.. Введение в метод конечных элементов: Пер. с англ. – М.: Мир, 1981. – 304с.
  5. Carlos A. Felippa. Introduction to Finite Element Methods. – Department of Aerospace Engineering Sciences and Center for Aerospace Structures University of Colorado, Boulder. – 2001.
  6. Kyran D. Mish, Leonard R. Herrmann, LaDawn Haws. Finite Element Procedures in Applied Mechanics (попадалось где-то в инете).
  7. Зенкевич О., Морган К. Конечные элементы и аппроксимация. – М.: Мир, 1986. – 318 с.
  8. Зенкевич О., Чанг И. Метод конечных элементов теории сооружений и в механике сплошных сред. – М.: Недра, 1974. – 240 с.
Ссылки:
  • http://www.fea.ru/ ...Сайт FEA.RU посвящен актуальным проблемам конечно-элементной механики и компьютерного инжиниринга (CAE), МКЭ и расчётам на прочность;
  • http://www.cae.ru/ форум CAD и CAE-систем в том числе теоретические и прикладные аспекты КЭ моделирования и решения задач механики деформируемого твердого тела. Механика конструкций, машин, сооружений и установок;
  • http://www.engr.usask.ca/~macphed/finite/fe_resources/fe_resources.html - Мощный каталог ресурсов, касающихся МКЭ;
  • http://www.isib.cnr.it/~secchi/EdMultifield/ - Сайт программы для конечноэлементных расчетов с неплохим описанием метода.

LinkWithin

Related Posts Plugin for WordPress, Blogger...