计算力学(五)有限元简介和插值函数
有限元法简介
里兹法的缺点
- 在求解区域形状不规则时,要在整个求解区域寻求一合适的容许函数是十分困难的
- 里兹法中的待求未知量是容许函数中包含的许多待定系数,而这些待定系数的物理意义并不明确
里兹法的缺点可通过分片插值函数来克服.即将整个求解区域分成若干个叫做单元的子区域,这些子区域的顶点叫做节点；并且假定单元和单元间只通过节点连结在一起.在每个单元上,容许函数都表示成单元节点物理量的插值函数的形式.由于划分的单元个数总是有限的,这种分片插值的里兹法叫做有限单元法(Finite Element Method),简称有限元(FEM)
有限元法的优点
- 对任意求解区域能较好的适应
- 采取分段函数作为容许函数,比较容易建立
- 近似函数中的系数有比较明确的物理意义
有限元法的一般步骤
- 离散并选择元素类型
- 位移函数
- 定义应变/位移和应力/应变关系
- 导出单元刚度矩阵和方程
- 将单元方程组合得到整体方程或总方程，并引入边界条件
- 求解未知自由度(或广义位移)
- 求解单元应变和应力
- 解读结果
插值函数
由于我们最终求得的是离散的网格点上的近似解,其他区域的点就需要依据求得的数据进行插值得到.
一维插值
一次插值
在一维区域中,假设已经求得在\(x=x_1,x=x_2\)处的函数值分别为\(y_1,y_2\),考虑最简单的情况,在\([x_1,x_2]\)区域上的函数值 \[ y(x)=f(y_1,y_2)=N_1y_1+N_2y_2\\ y(x_1)=y_1\quad, \quad y(x_2)=y_2\\ \] 不妨令\(\phi\)为最简单的多项式函数,\(x\)的最高次数最少1才能满足两个边界条件, \[ y(x)=f(y_1,y_2)=N_1y_1+N_2y_2=ax+b\\ y_1=ax_1+b\quad, \quad y_2=ax_2+b\\ \begin{align} y(x)&=\frac{y_2-y_1}{x_2-x_1}x+\frac{-x_1y_2+x_2y_1}{x_2-x_1}\\ &=\frac{x-x_2}{x_1-x_2}y_1+\frac{x-x_1}{x_2-x_1}y_2\\ \end{align}\\ \] 为了方便不同区间位置/大小的基函数有统一的形式,将区间转化到局部坐标系\([-1,1],\xi_1=-1,\xi_2=1\)得 \[ y(\xi)=\frac{1}{2}(1-\xi)y_1+\frac{1}{2}(1+\xi)y_2=\sum_{i=1}^{2}N_iy_i\\ \] ### 二次插值
有两种提高插值的精度的方法
- 划分更小的插值区间(h方法)
- 提高插值函数的阶次(p方法)
实践表明,\(p\)方法的收敛性大大优于\(h\)方法.在局部坐标系加入一个数据点\(y\bigg|_{\xi_3}=y_3\),其中\(\xi_3=0\)得到新的插值函数 \[ y(\xi)=\frac{1}{2}\xi(1-\xi)y_1+\frac{1}{2}\xi(1+\xi)y_2+(1-\xi^2)y_3=\sum_{i=1}^{3}N_iy_i\\ \] 也可以通过增加更多区间上的点来增加插值的精度,但二次插值是性价比较高的做法.上面解方程的步骤也可以用拉格朗日插值函数代替,只是我觉得有些复杂不太直观的情况还是解方程来的直接一些.
一次和二次插值的统一形式
\[ \begin{align} 若节点3存在&N_3=1-\xi^2\quad否则N_3=0\\ &N_1=0.5(1-\xi)-0.5N_3\\ &N_2=0.5(1+\xi)-0.5N_3\\ \end{align} \]
可以看出上面不管是一次插值还是二次插值都没考虑区间边界函数关系,因而在整体上只有函数值连续,我们可以要求边界处的导数也连续,得到更光滑的结果,这是Hermite插值,在梁、板、壳的有限元分析中有很大价值,具体操作大同小异就不写了.
插值函数的性质
- \(\displaystyle\sum_{i=1}^3N_i=1\)
- \(N_i(\xi_j)=\delta_{ij}\)
二维插值--矩形单元
双一次Lagrange插值函数
二维的局部坐标系变为了\(\xi\in[-1,1],\eta\in[-1,1]\)的正方形区域,
四节点插值的结果 \[ N_1(\xi,\eta)=0.25(1-\xi)(1-\eta)\qquad N_2(\xi,\eta)=0.25(1+\xi)(1-\eta)\\ N_3(\xi,\eta)=0.25(1+\xi)(1+\eta)\qquad N_4(\xi,\eta)=0.25(1-\xi)(1+\eta)\\ \] 可以看出二维拉格朗日插值基函数,可以由两个一维拉格朗日插值基函数相乘得到 \[ N_1(\xi,\eta)=\overline{N_1}(\xi)\overline{N_1}(\eta)\qquad N_2(\xi,\eta)=\overline{N_2}(\xi)\overline{N_1}(\eta)\\ N_3(\xi,\eta)=\overline{N_2}(\xi)\overline{N_2}(\eta)\qquad N_4(\xi,\eta)=\overline{N_1}(\xi)\overline{N_2}(\eta)\\ \] 其中\(\overline{N_i}(\xi)\)是一维一次拉格朗日插值基函数.4-9点的双一次Lagrange插值函数以此类推.最后得到\(u(\xi,\eta)=\displaystyle\sum N_i(\xi,\eta)u_i\)
Serendipity插值
类似一维一次和二次插值的统一形式,可以构造出8节点正方形二维的Serendipity插值 \[ 若节点5存在\quad N_5=0.5(1-\xi^2)(1-\eta)\quad 否则\quad N_5=0\\ 若节点6存在\quad N_6=0.5(1+\xi)(1-\eta^2)\quad 否则\quad N_6=0\\ 若节点7存在\quad N_7=0.5(1-\xi^2)(1+\eta)\quad 否则\quad N_7=0\\ 若节点8存在\quad N_8=0.5(1-\xi)(1-\eta^2)\quad 否则\quad N_8=0\\ N_1=\overline{N_1}-0.5N_8-0.5N_5\qquad N_2=\overline{N_2}-0.5N_5-0.5N_6\\ N_3=\overline{N_3}-0.5N_6-0.5N_7\qquad N_4=\overline{N_4}-0.5N_7-0.5N_8\\ \] 其中\(N_1\sim N_4\)是必须的,\(\overline N_i\)是是二维双一次Lagrange插值基函数.
插值函数的性质
- \(\displaystyle\sum_{i=1}^8N_i(\xi,\eta)=1\)
- \(N_i(\xi_j,\eta_j)=\delta_{ij}\)
二维插值--三角单元
对于三角形单元内点P的位置,可以用整体\(x,y\)坐标表示,也可以用面积坐标表示.下面定义面积坐标 \[ L_1=\frac{\Delta_1}{\Delta}\quad L_2=\frac{\Delta_2}{\Delta}\quad L_3=\frac{\Delta_3}{\Delta}\quad\\ \] 其中\(\Delta_i\)代表上图中各个三角形的面积,\(\Delta=\Delta_1+\Delta_2+\Delta_3\).显然\(L_1+L_2+L_3=1\),三个面积坐标中只有两个是独立的.下面直接给出6节点三角形插值函数 \[ \begin{align} 节点4存在\quad &N_4=4L_1L_2\quad 否则\quad N_4=0\\ 节点5存在\quad &N_5=4L_2L_3\quad 否则\quad N_5=0\\ 节点6存在\quad &N_6=4L_1L_3\quad 否则\quad N_5=0\\ &N_1=l_1-0.5N_4-0.5N_6\\ &N_2=l_2-0.5N_5-0.5N_4\\ &N_3=l_3-0.5N_6-0.5N_5\\ \end{align}\\ \]
三维插值
三维就不写了,在代码中给出20节点正六面体高阶单元.
补充
学了上面的之后,很自然的想到可以增加一个单元内节点的数量,从而得到更光滑/更精确(也许)的结果,通过待定系数解微分方程就能得到对应的插值函数,这里只是单纯的对称的增加边节点和内部节点数量,不涉及更高端的操作,据说上这门课的老师还发明了更高效/精确的单元,(我写这些只是之前自己研究过一点,不写出来感觉可惜了).但观察8节点正方形二次插值可知,插值函数中有\(\xi^2\eta,\xi\eta^2\)项,没有\(\xi^3,\eta^3\)项,同样的4节点中有\(\xi\eta\)项,没有\(\xi^2,\eta^2\)项,以此类推发现16点插值的时候最高次的项只有\(\xi^3\eta^3\), 而对应25点插值的时候最高次的项有\(\xi^4\eta^4\),对应就是下图项数三角形中梭形的部分
对应三维情况,也有类似的图
Periodic Table of the Finite Elements (umn.edu)列出了一些单元的形式,做的非常好
最后给出代码链接
#! https://zhuanlan.zhihu.com/p/693954230