c语言实现有限域模多项式_研一：王勖成-有限单元法第三、四章学习（二）

写在前面：本人小白，刚刚开始学习拓扑优化和C++相关知识，如果写的有错误请大家不吝赐教，也欢迎私信讨论~

这篇文章是学习王勖成的《有限单元法》后梳理的知识，所有截图来自该书，也可以理解为抄书hhh

一、弹性力学基本方程的矩阵形式

对于三维问题，弹性力学基本方程可以写成如下形式。

1、 平衡方程

弹性体再V域内任一点沿坐标轴x,y,z方向的平衡方程为

其中A 是为微分算子

2、 几何方程

几何方程的矩阵形式为

其中L为微分算子

3、 物理方程——应力应变关系

弹性力学中应力-应变之间的转换关系也称弹性关系，对于各项同性材料，应力可表示为

其中D称为弹性矩阵，它完全取决于弹性体材料的弹性模量E和泊松比v。

物理方程的另一种形式是

其中C是柔度矩阵，C=D-1，它和弹性矩阵是互逆关系。

4、 边界条件

包括力的边界条件和几何边界条件，此处暂不展开。

二、单元和插值函数的构造

着重系统地讨论了利用适合不同单元类型的局部坐标系建立各自单元插值函数的方法，为之后在规范域内采用标准化的步骤进行积分打下基础。

1、 一维单元

一维单元书上讲了拉格朗日单元和Hermite单元两种插值函数构造形式，我感觉拉格朗日单元重要一点。

对于具有n个结点的一维单元，单元内的长函数可以插值表示为

其中插值函数N 的性质不再赘述，关于此处插值函数N的构造，为了避免繁琐的推导，此处直接采用熟知的拉格朗日插值多项式，n个结点的一维单元就采用n-1次拉格朗日插值多项式：

如果放到自然坐标内，要求位移插值函数就要引入五晾干的局部坐标

其中l代表单元的长度，利用局部坐标可将插值函数构造为

这种无量纲表示为今后再自然坐标内积分打下基础。

2、 二维单元

对于二维的三角形单元，引入面积坐标为其自然坐标，下图为面积坐标内三角形单元（线性变化）的面积坐标图。

三角形内部或者边上的任意一点可由其面积坐标来表示即P（Li，Lj，Lm），其中

A是三角形的面积。由于其面积坐标与三角形的具体形状和其在总体坐标x，y中的位置无关，因此它是三角形的一种自然坐标。

对于二次变化的三角形单元，它拥有6个结点，各结点坐标和直线的方程标注如下

其各个结点对应的插值函数用划线法构造：

对于三次单元等高次单元，也可以使用划线法构造其插值函数。

通常情况下，二维的矩形单元比三角形单元更加方便有效。拉格朗日矩形单元的插值函数是利用两个坐标方向适当方次拉格朗日多项式的乘积得到的：

但是这种类型的单元存在一定缺点，因为其出现了随插值函数的方次增高而增加的内部结点，从而增加了单元的DOFs，而这些自由度的增加通常并不能提高单元的精度。这就引出了接下来更为重要的Serendipity四边形单元。

变结点数的Serendipity四边形单元可以实现同一单元的不同边界有不同数目的结点，这样可以进一步实现不同阶次单元之间的过渡，从而可能在求解的不同区域采用不同精度的单元。其插值函数的一般构造方法如下

3、 三维单元

与二维单元的面积坐标相对应的，三维单元的自然坐标是体积坐标，其构造方式与二维单元相对应。

4、 阶谱单元（hierachical element）

在自适应有限元分析中，当发现利用较低阶单元对同一问题进行分析的精度不满足给定的要求时，可能需要在单元网格划分不变的条件下提高单元的阶次。而阶谱单元可以实现以上目的，还可以让已经形成的低阶单元的刚度等特性矩阵保持不变地仍被利用，提高了运算效率。

比如在一维阶谱单元中，添加一个中间结点从而提高阶次使之成为二次单元时，将结点参数改写如下

这样以来，当单元从线性升阶为二阶单元时，原来的线性单元的相应部分可以保持不变地被继续使用，从而达到节省程序编写和运算时间的目的。但是在阶谱单元中，随着阶次的升高，其添加结点的结点参数的值将逐渐失去其物理意义，而且原来标准型C0型单元插值函数所具有的性质现在也不再保持了，但是单元内函数φ的插值表示的实质并未改变。

三、等参元和数值积分

上面给出了建立于局部的自然坐标内各类常用的单元及其插值函数的构造方法，在实际的工程问题和物理问题中，需要通过是适当的方法把上面说的规则形状的单元转化为其边界为曲线或者曲面的相应单元。此处最普遍采用的是等参变换，就是单元的几何形状和单元内的场函数采用相同树木的结点参数和相同的插值函数进行变换。

1、 等参变换

等参变换是为了将自然坐标中的几何形状规则的单元转换成总体（笛卡尔）总表中的几何形状扭曲的单元，前者成为母单元，后者称为子单元。

2、 单元矩阵的变换

有限元求解方程中需要用到单元体积内和面积内的积分，此处希望能在自然坐标内按照规格化的数值积分方法进行积分。

（1） 导数之间的变换

式中J称为雅克比矩阵，它可以显式地表示为自然坐标的函数。

（2） 体积微元、面积微元的变换

在笛卡尔坐标系内所形成的体积微元可以表示为：

通过以上坐标变换关系式以后，有限元问题中的积分式可以最终变换到自然坐标系的规则化域内进行：

3、 等参变换的条件和等参单元的收敛性

等参变换的条件是雅克比行列式|J|不得为零，进一步由雅克比行列式的公式推出：应防止任意的两个结点退化为一个结点，还应该防止元素的过分弯曲。

有限元分析中的收敛条件是元素是协调和完备的。为了保证协调，相邻单元的公共边或面上应该有完全相同的结点，同时它们采用相同的插值函数。为了满足完备性，要求插值函数中包含完全的线性项，即一次多项式。

4、 等参元用于分析弹性力学问题的一般格式

等参元通常以位移作为基本未知量，其一切计算都将在自然坐标系的母单元内进行，系统的方程需要做两方面更改：积分变量（取自然坐标）及积分限。

其中：

5、 数值积分方法

一维数值积分的基本思想是构造一个多项式去近似被积函数。所构造的多项式由积分点确定

此处介绍了两种积分：Newton-Cotes积分和高斯积分，其中高斯积分精度较高，应用广泛。在高斯积分的积分方案中积分点ξi不是等间距分布，其位置由下述方法确定：首先定义n次多项式P（ξ）

再由下列条件确定n个积分点的位置

被积函数F（ξ）可由2n-1次多项式φ（ξ）来近似，即

带入积分得到与Newton-Cotes积分相同形式的结果：

此处虽然形式上相同，但是实质上由区别：在高斯积分中φ（ξ）为2n-1次多项式，因此精度较高，又其积分点由方程确定，并非等间距的，优化了积分点的位置，所以精度高。

6、 等参元计算中数值积分阶次的选择