基于Mathematica的等参母单元形函数及其偏导数的推导

本期推文木木将带着大家用强大的符号计算软件Mathematica简单推一推有限元中的基本公式。

为什么用Mathematica推导？

为什么会有这种想法呢？木木最近在做UEL子程序开发的时候发现网上的子程序中，很多在求解母单元形函数偏导时直接把导数形式写上去，没有对其进行推导，当然，书中已然有很多等参元形函数及其偏导的形式，但是只限于经典的一些等参单元，比如说平面四节点、八节点、六面体八节点等一些常见的单元，那如果我们以后想开发自己的单元，节点形式在书中找不到，有时借助文献中的高阶单元，相信大多数情况下，文献中给出的高阶单元形函数是缩写指标形式，偏导更是需要自己推导，这时候就感受到推导公式的重要性了，以汪礼顺老先生在《工程力学》中发表的关于多结点等参单元形函数确定方法为例，

“  

多结点等参单元形态函数的确定——汪礼顺

如下图所示，为六面体32节点等参单元：

形函数的形式是这样的：

为此木木先抛砖引玉，利用当今符号计算强大的软件——Mathematica对平面八节点（Q8）单元以及六面体八节点（C3D8）单元的形函数进行展开、偏导数形式进行推导。

平面八节点（Q8）单元

可以先把该单元Plot一下，更加的直观。

Clear["Global`*"]
data = {{-1, -1}, {0, -1}, {1, -1}, {1, 0}, {1, 1}, {0, 1}, {-1, 
    1}, {-1, 0}, {-1, -0.99}};
Show[ListPlot[
  data -> {"节点1", "节点5", "节点2", "节点6", "节点3", "节点7", "节点4", "节点8", 
    "节点1"}, PlotStyle -> {Red}],
 ListLinePlot[data], AxesStyle -> Arrowheads[.02], 
 AxesLabel -> {"s", "t"}
 ]

形函数推导

接下来就进入正题了

由于代码片不支持上下标形式，所以会出现排版混乱的情况，木木就以截图的形式展现给大家。

输出结果：

大功告成，形函数展开形式可与上图对比，只与偏导形式，可以自行推导进行验证。

六面体（C3D8）单元

形函数推导

结果如下：