FEM07-从虚位移原理到有限元

3D结构问题

3D实体结构如图:

外载荷：

在XYZ坐标系中，相对于无载荷状态的位移U：

位移U产生的应变和应力：

需要分析的问题如下：

已知：实体的几何形状，外载荷，约束条件，材料的应力-应变关系和初始应力。

计算：位移、应变和应力。

目前只考虑线性问题，满足线性假设：

1、位移足够小，以至于可以基于无载荷状态建立平衡关系。

2、材料应力-应变关系可以随空间位置X、Y、Z变化，但在其他任何情况都保持不变。（如：C不随应力状态改变）

采用有限元方法，我们将实体离散为有限单元的组合。假设每个单元内的位移分布是所有N个单元节点位移的函数。

对于单元m，位移函数表示如下：

单元的位移函数u，可以在便于计算的单元局部坐标系xyz中定义（区别于大写字母表示的全局坐标系XYZ）。

H^(m)是位移插值矩阵（displacement interpolation matrix）。插值函数的图形很像一个帽子，也称为“Hat Function”。插值函数具备以下特征：H_i在节点i处取值为1，而在其他所有节点处都是0. （Takes a unit value at node i , and is zero at all other nodes.）

Û是所有节点的3个全局位移分量所构成的向量，它是一个3N维的向量。对于梁和板壳，还包括转角。虽然Û中包含了所有节点位移，但对于给定的单元，只有该单元的节点的位移会影响这个单元内的位移和应变分布。

推导单元应变：

B^(m)是应变-位移矩阵，由H^(m)的微分及其行的组合得到。

单元的应力由材料本构关系给出，并考虑初始应力：

在"FEM05-虚位移原理"中，介绍了虚位移原理的方程。现将其写成单元组合的形式（分别对每个单元运用虚位移原理，再求和）：

对于虚位移和虚应变，采用与真实的位移和应变相同的假设。代入上述方程，得到：

其中，K是总体刚度矩阵，它是由各个单元刚度矩阵K^(m)组装得到的。总体载荷向量R，包括体积载荷R_B、面载荷R_S、初始应力R_I的影响以及节点集中载荷R_C。需要注意：

1、所有要相加的矩阵的维数都是相同的。

2、单元自由度等于全局自由度。如果不一致，组装前需要进行转换。

例题

1、离散化，将结构分为两个单元。

2、单元的位移插值函数。

对于该问题所采用的单元，每个单元由两个节点（A和B）组成，长度为L，单元坐标系定义如下：坐标原点在A节点处，x方向由A指向B（与全局坐标系X方向一致）。只考虑轴向拉伸，所以只需要一个坐标方向。

假设单元内的位移u=H_A*U_A+H_B*U_B，且满足：

H_A在A节点处等于1，在其他节点处等于0;

H_B在B节点处等于1，在其他节点处等于0.

采用线性插值函数，可以取H_A=1-x/L和H_B=x/L.

3、分别写出各个单元的位移插值矩阵H^(m)、应变-位移矩阵B^(m)以及材料属性矩阵^(m)。

4、计算刚度矩阵K。

5、载荷向量R只包含节点集中载荷，R^T=[0,0,100]。

6、根据KU=R，求出各个节点的位移。

7、计算各个单元的应力。

上述采用线性插值的有限元方法得到的结果，与“FEM04-里玆法（Ritz Method）”中采用分段一次函数作为试函数的结果一致。主要区别在于：里兹法的试函数，虽然进行了分段，但每一段都是在全局坐标中定义的；而有限元中每个单元的插值函数，可以分别在单元局部坐标系中定义。

下面划分更多的单元，来求解这个例题。

将横截面变化的部分，分为4个单元。

分别写出各个单元的位移插值矩阵H^(m)、应变-位移矩阵B^(m)以及横截面积A^(m)。：

计算刚度矩阵K：

KU=R，其中U₁=0，通过矩阵运算得到位移和应力。（矩阵计算器：https://matrixcalc.org/en/）

对比不同单元数量的求解精度：

Reference:

• Finite Element Procedures by Klaus-Jürgen Bathe 轩建平 译
• A First Course in Finite Elements by Jacob Fish & Ted Belytschko