用CalculiX学习比赛

通过了折衷方案-您可以结婚!

引言


在复杂系统的设计中,有限元方法(FEM或FEM,他们已经在国外使用)牢固地进入了工程计算的实践。 在很大程度上,这与力学的强度计算有关。 该方法的应用(由相应的软件实现)大大减少了最终设备的开发周期,从而消除了在使用基于sopromat和结构力学方法的经典计算时进行必要的实验检查的麻烦。 迄今为止,已经开发了许多实现FEM的应用软件。 最前沿的是功能强大的ANSYS,在其两侧以及距离很远的地方-带有内置FEM模块的CAD系统(SolidWorks,Siemens NX,Creo Parametric,Compass 3D)。

CalculiX很强大,但是困难且难以理解。 我们会解决这个问题吗?



自然,FEM已渗透到教育领域-为了在实际任务中使用FEM,需要对相关专家进行培训。 在首都,在大型技术大学中,这方面的情况或多或少是正常的,在我们地区,例如,SFedU弹性理论系使用了相同的ANSYS。 但是在外围地区,在狭义的专业而不是富裕的大学中,情况令人遗憾。 而且很简单-一个工作场所的ANSYS成本约为200万卢布,而且需要多个地方。 不幸的是,并不是所有的大学都能负担得起30-40百万美元来组织一门计算机课程,以教授FEM的使用。

一种选择是在教育过程中使用免费软件。 幸运的是,可以使用这种软件。 但是,实际上没有使用俄语的资料。 为了纠正这种情况,我将把这篇文章专门介绍CalculiX的介绍。CalculiX是一个开放的免费软件包,旨在使用有限元方法解决固体可变形体力学以及流体和气体力学的线性和非线性三维问题。

1.什么是CalculiX,以及从何处获得。 Windows安装


CalculiX软件包是一组控制台实用程序,包括用于准备初始数据的预处理器,FEM求解器和用于处理结果的后处理器。 CalculiX既可以独立使用,也可以作为其他产品的一部分使用,其中FreeCAD势头强劲。 另一个问题是,CalculiX在我国仍然鲜为人知,有关此资源唯一一篇文章直接说明了这一点

我将特别介绍以下有关Windows工作的资料,这是包括教育机构在内最常用和最常用的资料。 而且,在其中使用许多免费程序是一个公开的痛苦。

如果您从CalculiX官方网站上获取Windows软件包,则完全不清楚下一步该怎么做。 结合英语文档,许多人结束了该产品,然后转化为关于无法使用该产品的有毒评论。 在某种程度上,这是正确的-进入门槛确实很高。 但是我们仍然会尝试。

这个奇迹对于Windows有许多非官方的相对较新手友好的构建,其中包括Windows的bConverged CalculiX 。 我们从此处下载分发工具包,将其解压缩,然后使用标准的“进一步,进一步...”方法进行安装。 因此,安装并不构成特别的谜,对于没有经验的用户来说,安装是很容易的。 作为主要的工作环境,该软件包使用SciTE文本编辑器,该编辑器集成了对CalculiX组件的调用以及命令交互输入的可能性,并且看起来像这样(可单击)。



2.梁弯曲问题及其使用Sopromat方法的解析解


举一个简单的学生问题-弯曲钢梁,将其一端捏住,将垂直力F施加到另一根上。



问题的参数如下:F = 10 kN; l = 1 m是光束的长度; h = 0.1 m和b = 0.05 m是横截面的尺寸。 为简单起见,我们将不考虑梁自身的重量,因为梁重量为39 kg的梁明显小于所施加的载荷。 我们在梁截面中找到最大法向应力,并计算由于弯曲变形而引起的梁挠度。

任何没有妥协的学生都可以轻松解决这一问题。 为了不让贵族们感到尴尬,我将决定的所有细节包裹在剧透中

通过以下方法解决问题
该问题是静态可确定的,并且可以简化为最简单的设计方案


在没有不适当的困难的情况下,我们从静力学方程中找到关系的反应

\开一个ñ X = 0 Y - F = 0MF\,l=0 endalign


从哪里来 M=F\,lY=F 。 弯矩图和单个弯矩图(应用Mohr积分所必需)被简单地构造并显示在图中。 弯曲梁的最大法向应力为

 sigmax= fracMz\,y maxIz


在哪里 y max=h/2=0.05 m是从截面的端点到梁的纵轴的最大距离; Iz -相对于弯曲轴的几何惯性矩,等于

Iz= fracb\,h312


通过对特定数据的简单计算,我们得出最大正常电压为  sigmax=119 兆帕

我们使用Mohr积分计算弯曲期间梁的最大挠度

 Deltaz= int limitsl0 fracMz\, overlineMzEIz\,dx= fracFE\,Iz int limitsl0x2\,dx= fracF\,l33\,E\,Iz


其中E = 200 GPa是钢的杨氏模量。 特定值的计算得出  Deltaz=3.97 cdot103 ,米

对于那些懒得看不起扰流板的人,我将立即给出问题的答案:电子束部分的最大正常电压  sigmax=119 MPa,最大挠度为3.97毫米。 给出这些数字是为了与随后在CalculiX中解决此问题的过程所提供的结果进行比较。

3.几何和计算网格的准备


首先,您需要将有关零件的几何数据输入CalculiX。 是的,可以从CAD导出几何图形,就像在同一ANSYS中所做的那样,但是我们将经历酷刑并手动输入几何图形。 从bConverged套件中打开SciTE编辑器,然后输入以下文本

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

使用名称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)。 稍后将在网格生成中派上用场。 现在用我们的耳朵假装

seta lines l l1 l2 l3 l4
swep lines sweeplines tra 0 0 0.1 10

一队

seta [ ] [ 1] ... [ N]

将多个对象组合成具有给定名称的集合。 实际上,这类似于对象分组。 下一条命令

swep [ ] [ ] [ ] [, ] []

移动选定的对象集以形成新的对象集。 可移动对象被复制。 在这种情况下,点的运动形成线,线的运动-曲面,表面的运动-连续体积。 在我们的例子中,我们将一组线沿Z轴移动了0.1米,而所得的线则分为10段。 我们按F10 ... er,这是什么?



嗯,黑屏...很容易修复,只需在脚本末尾添加行

plot pa all
plus la all

这些命令告诉您绘制所有点(pa)并将所有线条(la)添加到显示中,然后我们得到此结果



现在让我们基于我们创建的线组来创建曲面

seta surfaces s A001 A002 A003 A004

在脚本的最后添加这些表面的显示

plus sa all




现在,我们将执行另一个移动,即沿着Y轴移动0.05米,将由位移形成的所有直线展开5个分段。

swep surfaces swepsurface tra 0.0 0.05 0.0 5

振作精神



按住鼠标左键,然后删除点和线的显示,可以旋转生成的图像,我们将看到一些可理解的内容



是的... CalculiX与大众用户所熟悉的通常的视觉概念相去甚远,但是尽管如此,我们还是建立了光束的几何形状。

几何,几何,但是对于网格生成,我们将采取下一步行动-删除所有绘图和加号命令,并将几何生成代码包装在seto和setc命令中,如下所示

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

-根据名为beam的几何体生成由平行六面体(he8)组成的网格。 现在将生成的网格物体打印到文件中

send beam abq

-将网格输出以ABAQUS FEM软件包的格式输出到名为beam.msh的文件(有这样一个专有的FEM计算软件包,而CalculiX理解其格式)



因此,生成了网格,您可以查看beam.msh文件,并在其中看到类似的内容

*NODE, NSET=Nbeam
1,0.000000000000e+000,0.000000000000e+000,1.000000000000e-001
2,0.000000000000e+000,0.000000000000e+000,9.000000000000e-002
3,0.000000000000e+000,1.000000000000e-002,9.000000000000e-002
4,0.000000000000e+000,1.000000000000e-002,1.000000000000e-001
5,1.000000000000e-002,0.000000000000e+000,1.000000000000e-001
6,1.000000000000e-002,0.000000000000e+000,9.000000000000e-002
7,1.000000000000e-002,1.000000000000e-002,9.000000000000e-002
8,1.000000000000e-002,1.000000000000e-002,1.000000000000e-001
9,0.000000000000e+000,2.000000000000e-002,9.000000000000e-002
10,0.000000000000e+000,2.000000000000e-002,1.000000000000e-001
11,1.000000000000e-002,2.000000000000e-002,9.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

显然,这是网格元素的顶点及其坐标的列表,其次是面的列表。 为了使整个事情看起来更漂亮,我们使用了交互模式CalculiX。 为此,请保持图形窗口处于活动状态 ,然后依次输入以下命令

plot f beam

-显示几何的所有面

view edge off

-关闭边缘显示

view elem

-打开网格元素的显示。 我们通过按Enter键完成每个命令的输入,输入的命令将显示在右下方的SciTE窗口中,如下所示



是的,您不能称它为超级方便,但是尽管如此,我们仍然可以看到生成的网格物体的图像。



我注意到创建几何图形时创建的所有中间点都变成了网格节点。 因此,我们得到了一个六角形网格,尺寸为100 x 10 x 5个节点,元素边缘尺寸为10 mm。 我们创建的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



4.设定限制



应用有限元法的一个重要步骤是对结构点的移动设置限制,即考虑到对其施加的约束。 在我们的情况下,梁的一端被夹住,可以假设其一端是完全固定的。 我们需要告诉求解器FE网格的哪些节点是静止的。 当打开beam.fbd文件时,我们按F10键,等待带有梁图像的窗口出现



在交互方式下,输入命令
rot -y
plot n beam

第一小组部署模型,以便Y轴远离我们,第二小组-包括FE网格的绘制节点(n)。 移动模型(按住鼠标右键)并缩放图像(按住鼠标滚轮),我们得到了这张照片



现在我们需要选择所有要定义为固定的节点。 为此,我们再次使用交互式数据输入模式。 我们招募团队

qadd fixed

这将开始创建一组称为fixed的节点。 图形窗口上的光标切换到元素选择模式-以带有小方块的箭头形式显示。 这样放置光标



然后按r键。 然后像这样放光标



然后再次按r。 因此,我们形成了一个矩形的选择区域,其对角线由按r标记的光标位置设置。 我们用这个矩形选择我们需要位于梁末端的节点



按a,然后按n,突出显示标记的节点。 一条脚印将出现在控制台窗口中,其中包含选定节点的列表(图片可单击)



输入q退出选择模式并命令

plus n fixed g

以绿色(g)显示固定组的节点。 现在我们可以看到哪些节点将包含在固定条件中。



现在,我们必须卸载这些节点作为约束文件,然后将其馈送到求解器的输入中。 为此,请键入命令

send fixed abq spc 123

-以ABAQUS(abq)格式的约束文件的形式卸载一组固定节点,从而限制所有组节点在所有三个自由度(1-X轴,2-Y轴,3-Z轴)上的移动。 结果,形成了文件fixed_123.bou,内容如下

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


-实际上,这是所有节点的枚举以及对给定节点的移动进行限制的自由度数。

5.负载分配


固定梁之后,我们将尝试加载它。 再次打开面和元素的显示

plot f beam
view edge off
view elem

调整图像的方向,以便我们可以看到光束松动端的上部



让我们切换到对象选择模式

qadd load

将光标放在所需的脸部上,然后按f



该工作面以紫色突出显示,并在控制台窗口中显示一条说明,该说明已添加到载荷集中。

qadd load
2541 e:3873 s:6 n= 5298 5310 5312 5300

按a结束设定的形成,按q退出选择模式。 我们向选定的面施加压力,从而产生10,000 N的合力。很容易计算出选定面的面积为1 cm 2 ,这意味着所需压力为10 8 Pa。 使用命令设置此负载

send load abq pres 1e8

-以ABAQUS格式在load.dlo文件中显示负载。 该文件看起来像这样

** Pressure based on load
3873, P6, 100000000.000000

显示了网格元素的编号,其面以及该面上的压力值。 因此,可以认为初始数据的准备已经完成。

6.输入数据的描述和求解器的启动


现在应将所有这些数据(网格,限制和负载)拉到FEM求解器的输入中,为此,我们形成了这种输入文件

光束输入
*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.3
*SOLID SECTION,ELSET=Ebeam,MATERIAL=EL
*STEP
*STATIC
*DLOAD
*INCLUDE,INPUT=load.dlo
*NODE FILE
U
*EL FILE
S
*END STEP


我将更详细地解释什么是什么。 文件的第一部分

*HEADING
Model: CalculiX Beam Input File for Habrahabr article
*INCLUDE,INPUT=beam.msh

设置任务描述,并包括一个带有CE网格beam.msh的文件。 下一部分形成边界条件-我们在fixed_123.bou文件中定义的那些关系

*BOUNDARY
*INCLUDE,INPUT=fixed_123.bou


我们不应忘记设置为弹性的材料,确定其杨氏模量和泊松比。 我们取结构钢的平均值

*MATERIAL,NAME=EL
*ELASTIC
2e11,0.3
*SOLID SECTION,ELSET=Ebeam,MATERIAL=EL

最后一部分设置任务的类型-静态载荷的计算以及来自load.dlo文件的那些载荷

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


检查我们在SciTE中是否有一个带有beam.inp文件的选项卡后,按Ctrl + F10,从而启动求解器。 我们不知疲倦地告诉我们CalculiX已经在那里为我们计算了一些东西。 精疲力尽,以免弄乱我放在扰流板下方的文字

控制台输出光束问题的求解器
********************************************************** **********

CalculiX 2.10版,版权所有©1998-2015 Guido Dhondt
CalculiX绝对不提供保修。 这是免费的
软件,欢迎您在以下位置重新分发该软件
某些条件,请参阅gpl.htm

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

您正在使用2016年5月23日星期一13:24:06制作的可执行文件

以下数字是估计的上限

数量:

节点数:6666
元素:5000
一维元素:0
二维元素:0
每个元素的积分点:8
每个节点的自由度:3
每个元素的层数:1

分布的面部负荷:1
分布体积载荷:0
集中载荷:0
单点约束:198
多点约束:1
所有多点约束中的项:1
领带约束:0
受循环约束约束的从属节点:0
预紧约束中的从属节点:0

套:2
所有条款的总和:18332

材料:1
每种材料和温度的常数:2
每种材料的温度点:1
每种材料的塑料数据点:0

方向:0
振幅:2
所有振幅的数据点:2
打印请求:0
转换:0
财产卡:0

步骤1

选择静态分析

降低MPC的等级

确定矩阵的结构:
方程数
19800
非零下三角矩阵元素的数量
655236

最多使用1 cpu(s)进行应力计算。

使用高达1 cpu(s)的对称刚度/质量贡献。

使用对称线轴求解器分解方程组
线轴最多使用1 cpu。

最多使用1 cpu(s)进行应力计算。

工作完成

7.后处理和决策分析


求解器获得的结果需要后处理器进行处理。 要调用它,请按Shift + F10并获得带有光束图像的图形窗口。 单击该窗口的左侧,在带有光束图像的框架外,并获得菜单



我们有什么兴趣? 梁截面中的应力-选择数据集->应力。 菜单将消失,但我们再次调用它,然后选择“数据集”->“实体”->“ Mises”。 结果,von Mises等效应力模式打开。



所以-关键时刻! 梁截面的最大等效应力为117 MPa,与折衷结果略有不同。 但是! 为了解决sopromat问题,我们没有考虑弯曲和剪切过程中的切向应力,而是仅计算了弯曲产生的法向应力。 偏转将发生什么? 转到菜单:数据集-> DISP和数据集->实体-> D3



我们观察到最大位移对应于梁的加载端,等于3.96毫米! 使用Mohr积分宏伟且与我们的计算相关。

通过简单的操作( 可在此处阅读) ,还生成了光束变形的动画。



得出结论


“呃,伙计,请稍等,下一步是什么?!” 嗯,在提到FEM(尤其是CalculiX)时,人们无法在一整篇文章中提到各种各样的问题。 这篇文章实在是太多而且很无聊。 其目的是用一种可理解的语言来解释两件事:

  1. 开源未通过FEM分析软件
  2. 学习和使用此软件并不像乍看起来那样困难

足够一篇评论文章吗? 我也这么认为 在准备本文时,使用了以下来源:

  1. Calculix有限元分析梁 -用作所介绍所有材料的基础。 由于此处添加了作者获得的经验,并且所有代码均由他在撰写本文时编写,因此这不是翻译,而是俄语教程
  2. 官方CalculiX手册

Gitlab上提供了示例代码。

最后,我注意到-我并不坚强,我在大学里没有妥协。 过了一会儿,生活(和爱!)迫使我知道它的基础。 因此,文本中可能存在错误,我正在等待有关恶意评论的内容,并且我承诺将考虑所有评论。

感谢您的关注!

Source: https://habr.com/ru/post/zh-CN423359/


All Articles