Хабрахабр

Изучаем сопромат с CalculiX

Сдал сопромат — можно жениться!

Метод конечных элементов (МКЭ или FEM, у них за рубежом) прочно вошел в практику инженерных расчетов при проектировании сложных систем. В значительной степени это касается прочностных расчетов механики. Применения этого метода, реализуемого соответствующим программным обеспечением существенно сокращает цикл разработки конечного устройства, позволяя исключить массу экспериментальных проверок, необходимых при использования классических расчетов на основе методов сопромата и строительной механики. На текущий момент разработана масса прикладного ПО, реализующего МКЭ. Во главе угла стоит мощный ANSYS, по бокам от него и в почетном удалении — CAD-системы со встроенным FEM-модулем (SolidWorks, Siemens NX, Creo Parametric, Компас 3D).

Исправим это?
CalculiX силен, но труден и непонятен.

В столицах, в крупных технических вузах обстановка в этой области более-менее нормальная, да и у нас в регионе тот же ANSYS применяется, например, на кафедре теории упругости ЮФУ. Естественно, МКЭ проник и в сферу образования — чтобы использовать его в реальных задачах, нужна подготовка соответствующих специалистов. И всё просто — ANSYS стоит порядка 2 млн. Но по периферии, в узко специализированных и не богатых университетах ситуация плачевна. К сожалению не все вузы могут позволить себе выложить 30-40 миллионов на организацию компьютерного класса для обучения применению МКЭ. рублей за одно рабочее место, а место требуется не одно.

К счастью таковое ПО имеется. Одной из альтернатив может служить применение в учебном процессе свободного ПО. Исправляя эту ситуацию, данную статью я собираюсь посвятить в введению в CalculiX — открытый, свободный программный пакет, предназначенный для решения линейных и нелинейных трёхмерных задач механики твёрдого деформируемого тела и механики жидкости и газа с помощью метода конечных элементов.
Однако, русскоязычных материалов по его использованию практически не существует.

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

При этом использование многих свободных программ в ней представляет собой откровенную боль. Я специально буду излагать весть нижеследующий материал применительно к работе в ОС Windows, как самой распространенной и применяемой в том числе и учебных заведениях.

В совокупности с англоязычной документацией это, для многих ставит крест на данном продукте, а потом выливается в ядовитые комментарии по поводу невозможности его использования. Если брать пакет для Windows с официального сайта CalculiX, то становится совершенно непонятно, что с ним делать дальше. Но всё же мы попытаемся. И отчасти это справедливо — порог вхождения действительно высок.

Качаем дистрибутив отсюда, распаковываем, и устанавливаем его стандартным методом «далее, далее...». Существует ряд неофициальных относительно дружественных к новичку сборок сего чуда под Windows, среди них bConverged CalculiX for Windows. В качестве основной рабочей среды в данном пакете используется текстовый редактор SciTE, в который интегрированы вызовы компонентов CalculiX, а также представляется возможность интерактивного ввода команд и выглядит оно примерно так (кликабельно) Установка, таким образом, не представляет собой особого таинства и вполне доступна неискушенному пользователю.

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

Для простоты не будем учитывать собственный вес балки, так как он, при массе балки 39 кг существенно меньше чем приложенная нагрузка. Параметры задачи таковы: F = 10 кН; l = 1 м — длина балки; h = 0,1 м и b = 0,05 м — размеры поперечного сечения. Найдем максимальное нормальное напряжение в сечении балки а так же вычислим прогиб балки в следствии деформации изгиба.

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

Решение задачи методами сопромата

Для тех кому лень смотреть под спойлер, зразу приведу ответ задачи: максимальное нормальное напряжение в сечении балки $\sigma_x = 119$ МПа, а максимальный прогиб составляет 3,97 миллиметра. Эти цифры приводятся для последующего сравнения с тем, что даст нам процедура решения этой задачки в CalculiX.
Прежде всего необходимо ввести в CalculiX геометрические данные о рассматриваемой детали. Да, тут возможен экспорт геометрии из CAD-ов, как это делается в том же ANSYS, но мы поёдем путем мучений и введем геометрию вручную. Открываем редактор SciTE из комплекта bConverged и набираем там следующий текст

25 0 0
pnt p3 0. pnt p1 0 0 0
pnt p2 0. 75 0 0
pnt p5 1. 5 0 0
pnt p4 0. Мы увидим примерно следующее 0 0 0

Сохраняем файл с именем beam.fbd и жмем F10 для запуска препроцессинга.

Команда pnt создает точку в пространстве с заданными координатами и синтаксис её таков

pnt [имя точки] [x] [y] [z]

Теперь соединим эти точки линиями, добавив в файл следующий текст

line l1 p1 p2 25
line l2 p2 p3 25
line l3 p3 p4 25
line l4 p4 p5 25

получив после нажатия F10 следующую картинку

Команда

Это пригодится потом для генерации сетки. line [имя линии] [точка 1] [точка 2] [число сегментов линии]

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

1 10

Первая команда seta lines l l1 l2 l3 l4
swep lines sweeplines tra 0 0 0.

[объект N]

объединяет несколько объектов в набор с заданным именем. seta [имя набора] [объект 1] ... Следующая команда По сути это аналог группировки объектов.

Перемещаемые объекты копируются. swep [имя объекта или набора объектов] [имя нового набора объектов] [вид перемещения] [вектор, задающий перемещение] [параметр]

выполняет перемещение выбранного набора объектов с образованием нового набора. В нашем случае мы смещаем набор линий lines на вдоль оси Z на 0. При этом движение точек образует линии, движение линий — поверхности, движение поверхностей — сплошные объемы. Жмем F10… э, а что это? 1 метра, при этом, образующиеся линии делятся на 10 сегментов.

Гхм, пустой экран… Исправляется это легко, достаточно добавить в конец скрипта строчки

plot pa all
plus la all

Эти команды указывают рисовать все точки (pa) и добавить к отображению все линии (la), после чего мы получаем такой результат

Теперь создадим поверхности, на базе созданного нами набора линий

seta surfaces s A001 A002 A003 A004

добавив в самый конец скрипта отображение этих поверхностей

plus sa all

Теперь выполним ещё один сдвиг, теперь уже вдоль оси Y на 0,05 метров, развивая все образующиеся при смещении линии на 5 сегментов.

0 0. swep surfaces swepsurface tra 0. 0 5

Получим что-то в духе 05 0.

Полученное изображение можно повертеть, зажав левую кнопку мыши, и убрав отображение точек и линий, мы увидим нечто вразумительное

Дааа… CalculiX далек от привычных массовому пользователю визуальных концепций, но тем не менее мы построили геометрию нашей балки.

Геометрия, геометрией, но для генерации сетки сделаем следующий ход — уберем все команды plot и plus и обернем код генерации геометрии в команды seto и setc, вот так

25 0 0
pnt p3 0. seto beam
pnt p1 0 0 0
pnt p2 0. 75 0 0
pnt p5 1. 5 0 0
pnt p4 0. 1 10
seta surfaces s A001 A002 A003 A004
swep surfaces swepsurface tra 0 0. 0 0 0
line l1 p1 p2 25
line l2 p2 p3 25
line l3 p3 p4 25
line l4 p4 p5 25
seta lines l l1 l2 l3 l4
swep lines sweeplines tra 0 0 0. Теперь эту геометрическую группу можно пропустить через генерацию сетки, указав после всего вышеупомянутого кода команды 05 0 5
setc beam

Эта пара команд объединяет всю созданную геометрию в некий блок геометрии с именем beam.

Теперь выводим сгенерированную сетку в файл elty beam he8
mesh beam

— генерирует сетку, состоящую из параллелепипедов (he8) на основе геометрии с именем beam.

send beam abq

— вывод сетки в файл с именем beam.msh в формате FEM-пакета ABAQUS (есть и такой проприетарный пакет МКЭ-вычислений и CalculiX его формат понимает)

Таким образом сетка сгенерирована, можно глянуть в файл beam.msh и увидеть там что-то подобное этому

000000000000e+000,0. *NODE, NSET=Nbeam
1,0. 000000000000e-001
2,0. 000000000000e+000,1. 000000000000e+000,9. 000000000000e+000,0. 000000000000e+000,1. 000000000000e-002
3,0. 000000000000e-002
4,0. 000000000000e-002,9. 000000000000e-002,1. 000000000000e+000,1. 000000000000e-002,0. 000000000000e-001
5,1. 000000000000e-001
6,1. 000000000000e+000,1. 000000000000e+000,9. 000000000000e-002,0. 000000000000e-002,1. 000000000000e-002
7,1. 000000000000e-002
8,1. 000000000000e-002,9. 000000000000e-002,1. 000000000000e-002,1. 000000000000e+000,2. 000000000000e-001
9,0. 000000000000e-002
10,0. 000000000000e-002,9. 000000000000e-002,1. 000000000000e+000,2. 000000000000e-002,2. 000000000000e-001
11,1. 000000000000e-002
.
.
.
.
*ELEMENT, TYPE=C3D8, ELSET=Ebeam
1, 1, 2, 3, 4, 5, 6, 7, 8
2, 4, 3, 9, 10, 8, 7, 11, 12
3, 10, 9, 13, 14, 12, 11, 15, 16
4, 14, 13, 17, 18, 16, 15, 19, 20
5, 18, 17, 21, 22, 20, 19, 23, 24
6, 5, 6, 7, 8, 25, 26, 27, 28
7, 8, 7, 11, 12, 28, 27, 29, 30
8, 12, 11, 15, 16, 30, 29, 31, 32
9, 16, 15, 19, 20, 32, 31, 33, 34

По-видимому это список вершин элементов сетки с их координатами, за которым следует список граней. 000000000000e-002,9. Для этого, оставляя активным графическое окно, вводим последовательно следующие команды Чтобы всё это дело выглядело более красиво, воспользуемся интерактивным режимом CalculiX.

plot f beam

— отображаем все грани геометрии

view edge off

— отключаем отображение ребер

Ввод каждой команды завершаем нажатием Enter, введенные команды отображаются в окне SciTE справа внизу, вот так view elem

— включаем отображение элементов сетки.

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

Таким образом мы получили гексагональную сетку размером 100 х 10 х 5 узлов, с размером ребра элемента в 10 мм. Замечу, что все промежуточные точки, что были созданы при создании геометрии стали узлами сетки. Созданный нами файл beam.fbd описывает геометрию задачи и процесс создания сетки.

Полный текст файла beam.fbd

seto beam
pnt p1 0 0 0
pnt p2 0.25 0 0
pnt p3 0.5 0 0
pnt p4 0.75 0 0
pnt p5 1.0 0 0
line l1 p1 p2 25
line l2 p2 p3 25
line l3 p3 p4 25
line l4 p4 p5 25
seta lines l l1 l2 l3 l4
swep lines sweeplines tra 0 0 0.1 10
seta surfaces s A001 A002 A003 A004
swep surfaces swepsurface tra 0 0.05 0 5
setc beam

elty beam he8
mesh beam
send beam abq

В нашем случае один из концов балки защемлен, а значить можно считать, что один из её торцов полностью неподвижен. Важным этапом применения МКЭ является задание ограничений перемещения точек конструкции, то есть учет наложенных на неё связей. Жмем F10 при открытом файле beam.fbd, дожидаемся появления окна с изображением балки Нам нужно сообщить решателю, какие узлы КЭ-сетки являются неподвижными.

Двигая модель (зажав правую кнопку мыши) и масштабируя изображение (зажав колесо мыши) добиваемся такой картинки В интерактивном режиме вводим команды
rot -y
plot n beam

Первая команда разворачивает модель так, чтобы ось Y смотрела от нас, вторая — включает отрисовку узлов (n) КЭ-сетки.

Для этого опять используем интерактивный режим ввода данных. Теперь надо выбрать все узлы, которые мы хотим определить как неподвижные. Набираем команду

Курсор на графическом окне переключается в режим выбора элементов — отображается в виде стрелки с маленьким квадратиком. qadd fixed

что начинает создание набора узлов с именем fixed. А затем ставим курсор так

и снова жмем r. Ставим курсор так

и жмем клавишу r. Выделяем этим прямоугольником нужные нам узлы, лежащие на торце балки

жмем a а затем жмем n, выделяя отмеченные узлы. Таким образом, мы сформировали область выбора прямоугольной формы, диагональ которой задается отмеченными нажатием r положениями курсора. В консольном окне появится портянка со списком выделенных узлов (картинка кликабельна)

Вводим q для выхода из режима выбора и команду

Теперь мы можем увидеть, какие узлы у нас будут входить в условие закрепления
plus n fixed g

для отображения узлов группы fixed зеленым цветом (g).

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

В итоге образуется файл fixed_123.bou, следующего содержания send fixed abq spc 123

— выгрузить группу узлов fixed в виде файла ограничений в формате ABAQUS (abq), ограничив перемещение всех узлов группы по всем трем степеням свободы (1 — ось X, 2 — ось Y, 3 — ось Z).

** BOUNDARY based on fixed
1, 1, ,
2, 1, ,
3, 1, ,
4, 1, ,
9, 1, ,
10, 1, ,
13, 1, ,
14, 1, ,
17, 1, ,
.
.
.

— по сути, это перечисление всех узлов и номера степени свободы, по которой ограничивается перемещение данного узла.

Снова включим режим отображения граней и элементов После того, как мы закрепили нашу балку, попробуем её нагрузить.

plot f beam
view edge off
view elem

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

Перейдем в режим выделения объектов

qadd load

Наводим курсор на нужную грань и жмем f

К выделенной грани приложим давление, дающее равнодействующую силу величиной 10000 Н. Грань подсвечивается фиолетовым, а в консольном окне появляется описание того, что она добавлена в набор load
qadd load
2541 e:3873 s:6 n= 5298 5310 5312 5300

Жмем a для окончания формирования набора, жмем q для выхода из режима выбора. Задаем эту нагрузку командой Нетрудно посчитать, что площадь выделенной грани равна 1 см2, а значит искомое давление равно 108 Па.

Файл выглядит так send load abq pres 1e8

— выводит нагрузку в файл load.dlo в формате ABAQUS.

000000

Указан номер элемента сетки, его грань и величина давления на эту грань. ** Pressure based on load
3873, P6, 100000000. Таким образом, подготовку исходных данных можно считать завершенной.

Все те данные — сетка, ограничения и нагрузки, следует теперь запулить на вход решателя МКЭ, для чего мы формируем входной файл такого вида

3
*SOLID SECTION,ELSET=Ebeam,MATERIAL=EL
*STEP
*STATIC
*DLOAD
*INCLUDE,INPUT=load.dlo
*NODE FILE
U
*EL FILE
S
*END STEP
beam.inp
*HEADING
Model: CalculiX Beam Input File for Habrahabr article
*INCLUDE,INPUT=beam.msh
*BOUNDARY
*INCLUDE,INPUT=fixed_123.bou
*MATERIAL,NAME=EL
*ELASTIC
2e11,0.

Первая секция файла Поясню подробнее, что тут к чему.

Следующая секция формирует граничные условия — те связи, что мы определили в файле fixed_123.bou *HEADING
Model: CalculiX Beam Input File for Habrahabr article
*INCLUDE,INPUT=beam.msh

задает описание задачи и включает файл с КЭ-сеткой beam.msh.

*BOUNDARY
*INCLUDE,INPUT=fixed_123.bou

Значения берем средние для конструкционных сталей Не следует забывать и о материале, который мы задаем упругим, определяя его модуль Юнга и коэффициент Пуассона.

3
*SOLID SECTION,ELSET=Ebeam,MATERIAL=EL

Последняя секция задает тип задачи — расчет статического нагружения и те нагрузки, из файла load.dlo который при этом действуют *MATERIAL,NAME=EL
*ELASTIC
2e11,0.

*STEP
*STATIC
*DLOAD
*INCLUDE,INPUT=load.dlo
*NODE FILE
U
*EL FILE
S
*END STEP

Получаем выхлоп, говорящий нам о том, что CalculiX что-то там для нас посчитал. Проверив, что мы имеем в SciTE открытой вкладку с файлом beam.inp жмем Ctrl + F10, запуская тем самым решатель. Выхлоб, дабы не загромождать текст привожу под спойлером

Консольная выдача решателя для задачи beam

************************************************************

10, Copyright© 1998-2015 Guido Dhondt
CalculiX comes with ABSOLUTELY NO WARRANTY. CalculiX Version 2. This is free
software, and you are welcome to redistribute it under
certain conditions, see gpl.htm

************************************************************

You are using an executable made on Mon May 23 13:24:06 2016

The numbers below are estimated upper bounds

number of:

nodes: 6666
elements: 5000
one-dimensional elements: 0
two-dimensional elements: 0
integration points per element: 8
degrees of freedom per node: 3
layers per element: 1

distributed facial loads: 1
distributed volumetric loads: 0
concentrated loads: 0
single point constraints: 198
multiple point constraints: 1
terms in all multiple point constraints: 1
tie constraints: 0
dependent nodes tied by cyclic constraints: 0
dependent nodes in pre-tension constraints: 0

sets: 2
terms in all sets: 18332

materials: 1
constants per material and temperature: 2
temperature points per material: 1
plastic data points per material: 0

orientations: 0
amplitudes: 2
data points in all amplitudes: 2
print requests: 0
transformations: 0
property cards: 0

STEP 1

Static analysis was selected

Decascading the MPC's

Determining the structure of the matrix:
number of equations
19800
number of nonzero lower triangular matrix elements
655236

Using up to 1 cpu(s) for the stress calculation.

Using up to 1 cpu(s) for the symmetric stiffness/mass contributions.

Factoring the system of equations using the symmetric spooles solver
Using up to 1 cpu(s) for spooles.

Using up to 1 cpu(s) for the stress calculation.

Job finished

Для его вызова жмем Shift + F10 и получаем графическое окно с изображением балки. Полученные решателем результаты требует обработки постпроцессором. Кликаем в левой части этого окна, за рамкой с изображением балки и получаем меню

Напряжения в сечениях балки — выбираем Datasets -> STRESS. Что нас интересует? В итоге включается режим отображения эквивалентных напряжений по фон Мизесу Меню исчезнет, но мы снова вызываем его и выбираем Datasets -> Entities -> Mises.

Максимальное эквивалентное напряжение в сечении балки равно 117 МПа, что немного расходится с результатом из сопромата. Итак — момент истины! Решая задачу сопромата, мы не учитывали касательные напряжения при изгибе и сдвиге, а посчитали только нормальное напряжения от изгиба. Но! Идем в меню: Datasets -> DISP и Datasets -> Entities -> D3 Что будет с прогибом?

Великолепно и соотносится с нашим расчетом с использованием интеграла Мора. Наблюдаем, что максимальное перемещение соответствует нагруженному концу балки и равно оно 3,96 миллиметра!

Путем несложных манипуляций, о которых можно прочитать тут, генерируется и анимация деформаций балки

Хм, чувак не может вместить в одну статью всего того многообразия вопросов, что поднимаются при упоминании МКЭ в общем, и CalculiX в частности. «Э-э-э, чувак, постой, а что же дальше?!». И цель её — доходчивым языком объяснить две вещи: Статья и так вышла объемистой и довольно занудной.

  1. Open Source не прошел мимо ПО для МКЭ анализа
  2. Изучить и использовать данное ПО не так сложно, как может показаться на первый взгляд

Достаточно для обзорной статьи? Думаю да. При подготовке статьи использованы следующие источники:

  1. Calculix FEA Beam — послужило основой всего изложенного материала. Так как здесь добавлен опыт, полученный автором и весь код написан им в процессе написания статьи — это не перевод а именно тьюториал на русском языке
  2. Официальный мануал по CalculiX

Код примера доступен на Gitlab.

Чуть позже, познать его самые основы меня заставила жизнь (и любовь!). В заключении отмечу — я не прочнист, у меня не было сопромата в университете. Так что ляпы, возможно, присутствуют в тексте, о чем жду злобных комментов и даю обещание учесть все замечания.

Благодарю за внимание!

Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть