开源FEM程序Calculix导入喷气发动机模型
概述
帖子介绍了开源fem程序calculix的基本情况，以部分程序为例讲解calculix的组织方式。给出了calculix在windows平台的安装过程，最后以calculix程序中的经典算例：喷气式发动机模型为例，展示该程序的显示效果。该程序的所有源代码均可从网络获取，如果因为网络的原因获取不了，请联系我。
calculix介绍
calculix是一个开源的有限元分析（FEA）软件，主要用于结构力学、热分析和流体力学等领域的数值仿真。它支持线性和非线性分析，广泛应用于工程和科学研究中。 以下是calculix的一些关键特点：
有限元求解器：calculix提供了强大的有限元求解功能，能够分析结构、热和流体等多种问题。 多种求解功能：支持静力学、动力学、热分析以及大变形、非线性材料等复杂问题。 开源许可：calculix在 GNU GPL 许可下发布，用户可以自由地修改和分发软件。 前后处理工具：calculix 提供了与前处理器（如PrePoMax）和后处理器（如Paraview）兼容的功能，方便用户进行模型的创建和结果的可视化。 并行计算：支持并行计算，可以在多个处理器上运行，提升求解效率。 输入文件格式：用户通过文本文件定义模型、边界条件和材 料属性，格式与其他有限元软件（如ABAQUS）相似，这意味着从abaqus导出的inp文件可以直接在calculix计算，这是calculix设计者的聪明之处，允许该程序与abaqus互动。
calculix的应用领域非常广泛，包括汽车、航空航天、建筑、机械制造等行业，尤其适用于需要处理复杂力学行为和大规模计算问题的工程应用。
拿到程序之后，根目录如下图所示：
下面逐个介绍每个文件夹的作用。
首先是与程序执行相关的文件：bin、cxx_windows、config、src，这几个文件负责gui界面程序的调用和链接，我们不用对他们有过多的了解。
doc是一些文档，里面是一些html文件，可以直接用浏览器打开，如下图所示
examples里面是例子，后续可以在软件的界面中打开，值得注意的是，这些例子与abaqus使用相同的计算文件：inp文件。
img文件是一些软件的图标，我们不用过多关注。你把里面的图片换成狗头都没有问题，只要别改变里面图片的名称。
当第一次安装的时候，点击bat文件，程序会开始安装，安装完成之后再点击该文件，程序会启动，出现gui界面。
calculix源代码
首先向大家说明，该程序是开源的，任何人都能在网上获取和使用，但请遵守开源协议，如果大家因为网络原因（你懂得）无法获取，可以私信我，我提供源代码和可执行程序。
calculix主要采用fortran编写，大约有1189个fortran子程序，很显然，如果想全部读懂，几乎是不可能的，也没必要，我们只需要从里面吸取我们想要的养分即可。这些程序如下图所示
这只是其中一小部分fortran代码，可以发现，该程序并没有按照功能进行分组，不像opensees那样，这会给我们阅读程序带来困难，我们只能通过子程序的名称来大致判断功能，然后再进入到子程序内部查看注释，如beamgeneralsections.f子程序，通过名称可以发现这是可beam单元截面属性相关的子程序，然后我们再打开子程序，如下图所示
上面图片中箭头1指的地方是一些声明，有关于程序本身的开源协议等等。箭头2则是这个子程序的注释，说明了这个子程序的作用，我查看了其他的子程序，每一个子程序的这个位置都有一个注释说明子程序的作用，这应该是一个惯例，方便我们理解程序。这里面的注释说明了这个子程序是用来读取inp文件的界面属性的，我们往后看，继续验证我们的想法，下面我截取一部分程序
do i=2,n
if(textpart(i)(1:9).eq.'MATERIAL=') then
material=textpart(i)(10:89)
elseif(textpart(i)(1:12).eq.'ORIENTATION=') then
orientation=textpart(i)(13:92)
elseif(textpart(i)(1:6).eq.'ELSET=') then
elset=textpart(i)(7:86)
elset(21:21)=' '
ipos=index(elset,' ')
elset(ipos:ipos)='E'
elseif(textpart(i)(1:8).eq.'SECTION=') then
if(textpart(i)(9:12).eq.'PIPE') then
section='PIPE'
ndprop=2
elseif(textpart(i)(9:11).eq.'BOX') then
section='BOX'
ndprop=6
elseif(textpart(i)(9:15).eq.'GENERAL') then
section='GENE'
ndprop=5
else
write(*,*)
& '*ERROR reading *BEAM SECTION: unknown section'
ier=1
return
endif
elseif(textpart(i)(1:8).eq.'OFFSET1=') then
read(textpart(i)(9:28),'(f20.0)',iostat=istat) offset1
if(istat.gt.0) then
call inputerror(inpc,ipoinpc,iline,
& "*BEAM SECTION%",ier)
return
endif
elseif(textpart(i)(1:8).eq.'OFFSET2=') then
read(textpart(i)(9:28),'(f20.0)',iostat=istat) offset2
if(istat.gt.0) then
call inputerror(inpc,ipoinpc,iline,
& "*BEAM SECTION%",ier)
return
endif
else
write(*,*) '*WARNING reading *BEAM SECTION:'
write(*,*) ' parameter not recognized:'
write(*,*) ' ',
& textpart(i)(1:index(textpart(i),' ')-1)
call inputwarning(inpc,ipoinpc,iline,
& "*BEAM SECTION%")
endif
enddo
熟悉abaqus的计算inp文件的家人们都能看出来，上面的程序涉及到“SECTION、MATERIAL”等等，这就是inp文件中声明和赋予单元材料属性的关键字。事实上，经过我后续学习，这确实是一个提取inp文件中材料属性的子程序。
上面展示了一个提取inp文件材料属性的子程序，下面我在给出一些 关于积分点的子程序供大家学习。
!
! contains Gauss point information
!
! gauss1d1: lin, 1-point integration (1 integration point)
! gauss1d2: lin, 2-point integration (2 integration points)
! gauss1d3: lin, 3-point integration (3 integration points)
! gauss2d1: quad, 1-point integration (1 integration point)
! gauss2d2: quad, 2-point integration (4 integration points)
! gauss2d3: quad, 3-point integration (9 integration points)
! gauss2d4: tri, 1 integration point
! gauss2d5: tri, 3 integration points
! gauss2d6: tri, 7 integration points
! gauss3d1: hex, 1-point integration (1 integration point)
! gauss3d2: hex, 2-point integration (8 integration points)
! gauss3d3: hex, 3-point integration (27 integration points)
! gauss3d4: tet, 1 integration point
! gauss3d5: tet, 4 integration points
! gauss3d6: tet, 15 integration points
! gauss3d7: wedge, 2 integration points
! gauss3d8: wedge, 9 integration points
! gauss3d9: wedge, 18 integration points
! gauss3d10: wedge, 6 integration points
! gauss3d11: wedge, 1 integration points
! gauss3d12: hex, 14 integration points (for c3d27)
! gauss3d13: hex, 2x5x5=50 integration points (for beams)
! gauss3d14: wedge, 1 integration point
!
! weight2d1,... contains the weights
!
!
real*8 gauss1d1(1,1),gauss1d2(1,2),gauss1d3(1,3),
& gauss2d1(2,1),gauss2d2(2,4),gauss2d3(2,9),gauss2d4(2,1),
& gauss2d5(2,3),gauss3d1(3,1),gauss3d2(3,8),gauss3d3(3,27),
& gauss3d4(3,1),gauss3d5(3,4),gauss3d6(3,15),gauss3d7(3,2),
& gauss3d8(3,9),gauss3d9(3,18),gauss3d10(3,6),gauss3d11(3,1),
& gauss2d6(2,7),gauss3d12(3,14),gauss3d13(3,50),gauss3d14(3,1),
& weight1d2(2),weight1d3(3),weight2d1(1),weight2d2(4),
& weight2d3(9),weight2d4(1),weight2d5(3),weight3d1(1),
& weight3d2(8),weight3d3(27),weight3d4(1),weight3d5(4),
& weight3d6(15),weight3d7(2),weight3d8(9),weight3d9(18),
& weight3d10(6),weight3d11(1),weight2d6(7),weight3d12(14),
& weight3d13(50),weight3d14(1),weight1d1(1)
!
gauss1d1=reshape(( /0.d0/),(/1,1/))
!
gauss1d2=reshape(( /
& -0.577350269189626d0,0.577350269189626d0/),(/1,2/))
!
gauss1d3=reshape(( /
& -0.774596669241483d0,0.d0,0.774596669241483d0/),(/1,3/))
!
gauss2d1=reshape((/0.d0,0.d0/),(/2,1/))
!
! the order of the Gauss points in gauss2d2 is important
! and should not be changed (used to accelerate the code
! for CAX8R axisymmetric elements in e_c3d_th.f)
!
gauss2d2=reshape((/
& -0.577350269189626d0,-0.577350269189626d0,
& 0.577350269189626d0,-0.577350269189626d0,
& -0.577350269189626d0,0.577350269189626d0,
& 0.577350269189626d0,0.577350269189626d0/),(/2,4/))
!
gauss2d3=reshape((/
& -0.774596669241483d0,-0.774596669241483d0,
& -0.d0,-0.774596669241483d0,
& 0.774596669241483d0,-0.774596669241483d0,
& -0.774596669241483d0,0.d0,
& -0.d0,0.d0,
& 0.774596669241483d0,0.d0,
& -0.774596669241483d0,0.774596669241483d0,
& -0.d0,0.774596669241483d0,
& 0.774596669241483d0,0.774596669241483d0/),(/2,9/))
!
gauss2d4=reshape((/0.333333333333333d0,0.333333333333333d0/),
& (/2,1/))
!
gauss2d5=reshape((/
& 0.166666666666667d0,0.166666666666667d0,
& 0.666666666666667d0,0.166666666666667d0,
& 0.166666666666667d0,0.666666666666667d0/),(/2,3/))
!
gauss2d6=reshape((/
& 0.333333333333333d0,0.333333333333333d0,
& 0.797426985353087d0,0.101286507323456d0,
& 0.101286507323456d0,0.797426985353087d0,
& 0.101286507323456d0,0.101286507323456d0,
& 0.470142064105115d0,0.059715871789770d0,
& 0.059715871789770d0,0.470142064105115d0,
& 0.470142064105115d0,0.470142064105115d0/),(/2,7/))
!......
上面的子程序是有关于数值积分程序的一部分，全部的程序有四百多行，我就不全放出来了，大家可以在calculix源代码中找到。
calculix下载与安装
首先进入链接：https://github.com/calculix/cae，然后在下图中点击downloads。
进入如下界面
这四个文件中，从上到下，第一个是linux系统的可执行文件，第二个是windows系统中的可执行文件，后面两个是源代码，里面的内容是一样的，只不过压缩方式不同，这里以windows系统为例讲解安装过程，下载第二个即可，如果你需要学习源代码，后两个文件随便选择一个解压即可。
解压之后，文件夹中会有“cae.bat”文件，点击这个文件，程序会进入安装状态，会安装一些工作依赖库，如：pyQt、paraview等等，我们不用干涉，程序安装结束后，calculix会自动打开，并显示默认的文件，如下图所示
上面有是三个窗口，最左边1号窗口是命令行，里面会随着你的操作刷新一些命令，这个窗口不用管，中间的2号窗口是程序的主界面，我们主要的操作就是在这个界面中，右边的3号窗口是显示窗口，里面实时显示我们的操作结果。
喷气发动机导入
每次打开软件，calculix都会导入默认的模型，那如何导入我们自己的模型呢？
如下图所示
在菜单栏依次点击1（file）和2（import model），软件会弹出如下界面
上面的1中显示的是默认工作路径中的模型文件，2中是我们选择的文件类型，我们要导入的喷气发动机是一系列零件组装的，因此这里选择组装体，如下图所示
导入需要花一些时间，导入的过程中，命令行会持续刷新一些命令，这里随便截取最后刷新的命令如下图
可以看到，在导入模型的最后时刻，calculix尝试修复一些病态单元，最后更新了一些东西，当显示“done”字样的时候就说明，模型导入成功。
下面放几张喷气发动机模型的图片供大家欣赏
上面就是喷气发动机的显示效果，还是比较震撼的！
总结
开源fem程序calculix采用fortran程序编写，该程序采用与大型商业软件abaqus相同的计算文件，即inp文件，这极大的方便了两者之间的交互，使得该软件可以无缝调用abaqus导出的计算文件。开发者提供了可执行程序和全部的源代码，任何使用者都可以从网络上获取，里面所有的源代码都可以查看，非常方便学习。
