[原创]：讨论怎样使有限元程序能自定义单元？

chenkong

在此讨论怎样能使有限元程序能自定义单元。我有一个想法，拿出来跟各位讨论，希望能听听大家的建议。附件包含一个word文档。

面向对象编程方法出现很多年了，一直没能给有限元计算程序带来多大变化。我想应该是和旧的fortran时代有限元程序说再见的时候了。
使用面向对象编程，程序将对所有单元都统一接口。再进一步，程序能允许用户自定义刚度矩阵。我的办法是，程序可以采用读脚本（比如文本文件，当然现在xml文件格式炒得挺热）的方式自定义新的单元。
比如用户定义一个刚度矩阵：

如上图所示的一个平面杆单元的刚度矩阵。程序能读入这个矩阵，格式可能会是：
{
{E*A/l, 0, 0, -E*A/l , 0, 0},
{0, 12* E*A/l^3, -6* E*A/l^2, 0, -12* E*A/l^3, -6* E*A/l^2},
{0, -6* E*A/l^2, 4* E*A/l, 0, 6* E*A/l^2, 2* E*A/l},
{-E*A/l, 0, 0, E*A/l , 0, 0},
{0, -12* E*A/l^3, 6* E*A/l^2, 0, 12* E*A/l^3, 6* E*A/l^2},
{0, -6* E*A/l^2, 2* E*A/l, 0, 6* E*A/l^2, 4* E*A/l }
}
这样，假如某个用户经过研究，认为矩阵第1行第1列的元素是1.1*E*A/l更符合他的要求，他只需要新建一个符合预定格式的文本文件，把E*A/l改成1.1*E*A/l就行了。这样做无疑使程序更能适应变化。比如说有的用户希望考虑剪切对杆单元的影响，这样他就可以自己定义一个6*6的矩阵，只要格式符合要求，程序就能读入这个新的单元。
这样做会有一些难度，比如3维杆单元应该是9*9的矩阵，怎样和平面杆单元兼容？直观的办法，我们可以定义矩阵元素代表的含义。比如，平面杆单元的矩阵项代表：
node1X, node1Y, node1θ, node2X, node2Y, node2θ
3维杆单元的矩阵项代表：
node1X, node1Y, node1Z, node1θx, node1θy, node1θz,
node2X, node2Y, node2Z, node2θx, node2θy, node2θz
这样，我们在3维杆单元刚度矩阵中表示平面杆单元的话就只需要在node1Z,node1θx, node1θy, node2Z,node2θx, node2θy项填0就可以了。
这样做也会带来新的难点。稀疏矩阵对角线上如果有0元素就会引起奇异，所以引入0元素得考虑处理的办法。我能想到的是对总刚矩阵做一次操作，删除0元素。
另外，很容易想到，3维三角形壳单元将会是18*18的刚度矩阵，3维四边形壳单元将会是24*24的刚度矩阵。所以，程序还得能够正确识别单元对应的各个位移项——也就是我们预先定义的矩阵元素含义。
欢迎各位能对我的想法提出更多的建议。


陈睿
2006-7-22