瞬态热传导有限元求解器开发

进展自研的工程壳有限元软件在有一搭没一搭的开发中，第一次跑通了前后处理。首先说说我搞这个软件的目的。实际上，杆、梁、板、膜、壳、实体单元求解器我都写过。已经在实际工程应用的是梁单元，我用它解决了一些叶轮机械设计的问题，集成为了一些专用软件。对于我心心念念的复合材料力学问题，用实体单元开发会快一些，但是考虑到目前存在的UMAT/VUMAT这些子程序与ABAQUS 铺层模块互斥的问题，我想搞一个自己的复合材料壳，支持批量定义铺层，还能够把我之前写的一些损伤、失效、疲劳的本构放进去，计算效率还远远高于实体单元，这样就具有了实际工程价值。所以我称它是工程壳有限元。为什么说用实体单元开发会快一些呢?在仿真中，壳的计算效率高，网格划分起来也简单，给人的感受是好像壳比实体要简单。然而在开发上，恰恰相反。实体单元的刚矩阵简单，而壳是由膜单元和板单元两部分合并构成的，尤其是三角形单元，其形函数非常复杂。看了不少论文专门做这个推导，有些甚至推导到最后还是错的。所以为了搞壳单元，被迫先去搞了膜和板单元。开发思路在搞求解器的同时，我的主要精力放在了视口 交互上。几年前我就会基于现有的库做三维显示了，如果只是显示一下模型和结果，这个对我不难。视口 交互但是既然要做工程软件，要有些实用性。我们在用其他商用软件的时候有个明显的体会，如果不能总视口中框选、点选、按角度选进行节点或者单元的拾取，我定义边界条件，或者修改其他设置的时候会有多么麻烦。由于节点的拾取频率最高，所以框选、点选、按角度选，都要做，除了能选，还要选中后高亮，还要可以选中后撤销指定的点。用惯了别的商用软件，感觉这些功能很基础，真到自己开发的时候，这些反而是最难的。为什么难呢？或者说开发中大型工业软件的难点在哪呢？有人说求解器，有人会说算法。在我看来，最难的是数据结构。比如一个单元在视口中是怎么表达的，在求解器输入中如何表达，在结果中如何表达，他们之间又是如何映射的。当你选中一个节点的时候，它怎么知道自己被选中了，它的编号是多少，你拿到这个数据存在内存，当用户二次查看的时候，你又如何让它二次高亮上去？所以，要弄清楚，这些模型的数据在不同的模块应该如何存储、如何转换、如何调用。定义set框选选中后高亮特征标记用户在定义边界条件以后。需要针对固支、位移、载荷、弯矩等边界条件，在指定节点位置进行标记显示。有一个非常有意思的现象，当我们在ABAQUS定义Fx、Fy等点载荷的时候，他会同时显示两个方向的载荷箭头，而不是显示合力方向。以前我不理解，现在我来理解了，我们看到的箭头应该不是通过向量绘制的。应该是把其内置的图形在节点位置进行显示，同时显示多个方向的载荷箭头，在代码中实现最简单，如果搞成合力方向，还要去计算偏转角度，显示的效果也不好。当然这是我猜的，反正我也采用这个思路。包含了载荷、位移、固支约束的标记同样的问题，当更改边界条件的时候，这些标记还要自动更新，二次打开工程的时候，还要能重建这些特征。还是数据结构的难题。这也是我花了很多时间的在视口上的原因。既然迟早要做，不如趁现在，再难也要做出来。特点我借鉴了ABAQUS set的思路，凡是要事先定义边界条件的地方，必须要要先定义成set。后面定义边界条件的时候，这些set会自动出现在下拉列表供选取，并且也通过“显示”按钮进行查看。求解器难点前面讲了交互的难，实际上求解器也是个难点，难在应力求解上。前面我们提到了，壳是膜单元和板单元合并成的。在求应力应变的时候，还需要根据位移，分布在膜和板里面求应力，最后再拼起来。和实体单元比，壳和梁单元都有一个问题，就是局部坐标系和整体坐标系的转换。在做刚度矩阵的时候，需要从局部到整体。在计算应变的时候，又要先把整体坐标系的位移转换到局部，然后在局部求应变。总之弯弯绕绕很多。仿真效果平面模型我们软件的位移结果某商用软件位移结果我们软件的应力结果某商用软件应力结果曲面叶片模型 我们软件的位移结果某商用软件位移结果我们软件的应力结果某商用软件应力结果来源：静界有限元