当前位置：首页 > news >正文

【小呆的力学笔记】非线性有限元的初步认识【三】

news 文章来源：https://blog.csdn.net/u012915522/article/details/129812872 2025/5/14 9:06:45

文章目录

- - - 1.2.2 基于最小势能原理的线性有限元一般格式
    - - 1.2.2.1 离散化
      - 1.2.2.2 位移插值
      - 1.2.2.3 单元应变
      - 1.2.2.4 单元应力
      - 1.2.2.5 单元刚度矩阵
      - 1.2.2.6 整体刚度矩阵
      - 1.2.2.7 处理约束
      - 1.2.2.8 求解节点载荷列阵
      - 1.2.2.9 求解位移列阵
      - 1.2.2.10 计算应力矩阵等

1.2.2 基于最小势能原理的线性有限元一般格式

现在市场上的有限元分析软件普遍都是基于位移法进行求解的，其数学原理就是上节的最小势能原理变分提法。回顾上节推导过程，首先需要根据位移条件确定可能位移的范围，其次根据假设的可能位移代入几何方程和本构方程，得到基于位移的应力应变的表达式，最后代入 $\delta II=0$ 的方程求解相关待定参数。线性有限元求解的一般格式基本上也是一样的，除此之外线性有限元会多一些其他步骤，比如结构的离散化、位移模式假定等。下图为线性有限元求解的一般格式。
（加图）

1.2.2.1 离散化

将结构进行离散化建模，是有限元分析工作中的第一步，某典型的结构如下图所示，其中结构为一个圆盘，将其离散成若干四边形网格单元，如下图所示。结构离散化合适与否一定程度上决定使用有限元分析方法解决工程问题的精度。
在这里插入图片描述

1.2.2.2 位移插值

（a）位移模式
在最小势能原理中首先需要根据位移条件确定可能位移的范围，但是在通常的分析中，我们为了方便处理，往往将可能位移用多项式来逼近，这个过程就是单元的位移模式选择。
我们以一阶二维三角形单元为例，二维三角形单元有6各位移自由度。用多项式拟合可能位移，如下式所示

$\begin{aligned} u(x,y)=\overline a_{0}+\overline a_{1}x+\overline a_{2}y\\ v(x,y)=\overline b_{0}+\overline b_{1}x+\overline b_{2}y \end{aligned}\tag{1-54}$

在这里插入图片描述

（b）形函数
多项式位移模式有六个待定参数，二维三角形单元有六个节点位移变量，所以六个待定参数可以表示成六个节点位移变量的形式，如下式，待定系数满足下式。
$\begin{aligned} \begin{bmatrix} u_1\\ u_2\\ u_3\\ \end{bmatrix}= \begin{bmatrix} 1 & x_1 & y_1\\ 1 & x_2 & y_2\\ 1 & x_3 & y_3\\ \end{bmatrix} \begin{bmatrix} \overline a_0\\ \overline a_1\\ \overline a_2\\ \end{bmatrix}, \begin{bmatrix} v_1\\ v_2\\ v_3\\ \end{bmatrix}= \begin{bmatrix} 1 & x_1 & y_1\\ 1 & x_2 & y_2\\ 1 & x_3 & y_3\\ \end{bmatrix} \begin{bmatrix} \overline b_0\\ \overline b_1\\ \overline b_2\\ \end{bmatrix} \end{aligned}\tag{1-55}$
那么系数为
$\begin{aligned} \begin{bmatrix} \overline a_0\\ \overline a_1\\ \overline a_2\\ \end{bmatrix}= \begin{bmatrix} 1 & x_1 & y_1\\ 1 & x_2 & y_2\\ 1 & x_3 & y_3\\ \end{bmatrix}^{-1} \begin{bmatrix} u_1\\u_2\\u_3\\ \end{bmatrix}, \begin{bmatrix} \overline b_0\\ \overline b_1\\ \overline b_2\\ \end{bmatrix}= \begin{bmatrix} 1 & x_1 & y_1\\ 1 & x_2 & y_2\\ 1 & x_3 & y_3\\ \end{bmatrix}^{-1} \begin{bmatrix} v_1\\ v_2\\ v_3\\ \end{bmatrix} \end{aligned}\tag{1-56}$
其中逆矩阵的求解如下
$\begin{bmatrix} 1 & x_1 & y_1\\ 1 & x_2 & y_2\\ 1 & x_3 & y_3\\ \end{bmatrix}^{-1}=\frac{[M_{ij}\cdot (-1)^{i+j}]^T}{A}\tag{1-57}$
其中 $M_{ij}$ 为原矩阵的余子式，典型如下，其乘上 $1)^{i+j}$ 称为代数余子式，矩阵的逆是其代数余子式组成的伴随矩阵除以该矩阵的行列式。

在这里插入图片描述
经过求解，系数矩阵的逆矩阵如下式
$\begin{aligned} \begin{bmatrix} 1 & x_1 & y_1\\ 1 & x_2 & y_2\\ 1 & x_3 & y_3\\ \end{bmatrix}^{-1}&= \begin{vmatrix} 1 & x_1 & y_1\\ 1 & x_2 & y_2\\ 1 & x_3 & y_3\\ \end{vmatrix}^{-1}\cdot \begin{bmatrix} \begin{vmatrix} x_2 & y_2\\ x_3 & y_3\\ \end{vmatrix} & -\begin{vmatrix} x_1 & y_1\\ x_3 & y_3\\ \end{vmatrix} & \begin{vmatrix} x_1 & y_1\\ x_2 & y_2\\ \end{vmatrix}\\ & & \\ -\begin{vmatrix} 1 & y_2\\ 1 & y_3\\ \end{vmatrix} & \begin{vmatrix} 1 & y_1\\ 1 & y_3\\ \end{vmatrix} & -\begin{vmatrix} 1 & y_1\\ 1 & y_2\\ \end{vmatrix}\\& & \\ \begin{vmatrix} 1 & x_2\\ 1 & x_3\\ \end{vmatrix} & -\begin{vmatrix} 1 & x_1\\ 1 & x_3\\ \end{vmatrix} & \begin{vmatrix} 1 & x_1\\ 1 & x_2\\ \end{vmatrix}\\ \end{bmatrix} \\=\frac{1}{A}\begin{bmatrix} a_1 & a_2 & a_3\\ b_1 & b_2 & b_3\\ c_1 & c_2 & c_3\\ \end{bmatrix} \end{aligned}\tag{1-58}$
那么，系数可以表达成
$\begin{aligned} \overline a_0=\frac{1}{A}(a_1u_1+a_2u_2+a_3u_3)\\ \overline a_1=\frac{1}{A}(b_1u_1+b_2u_2+b_3u_3)\\ \overline a_2=\frac{1}{A}(c_1u_1+c_2u_2+c_3u_3) \end{aligned}\tag{1-59}$
上式代入位移模式，可得
$\begin{aligned} u(x,y)&=\overline a_{0}+\overline a_{1}x+\overline a_{2}y\\ &=\frac{1}{A}(a_1u_1+a_2u_2+a_3u_3)+\frac{1}{A}(b_1u_1+b_2u_2+b_3u_3)x+\frac{1}{A}(c_1u_1+c_2u_2+c_3u_3)y\\ &=\frac{1}{A}(a_1+b_1x+c_1y)u_1+\frac{1}{A}(a_2+b_2x+c_2y)u_2+\frac{1}{A}(a_3+b_3x+c_3y)u_3\\ &=N_1(x,y)u_1+N_2(x,y)u_2+N_3(x,y)u_3 \end{aligned}\tag{1-60}$
同理，可得
$\begin{aligned} v(x,y)&=\overline a_{0}+\overline a_{1}x+\overline a_{2}y\\ &=\frac{1}{A}(a_1v_1+a_2v_2+a_3v_3)+\frac{1}{A}(b_1v_1+b_2v_2+b_3v_3)x+\frac{1}{A}(c_1v_1+c_2v_2+c_3v_3)y\\ &=\frac{1}{A}(a_1+b_1x+c_1y)v_1+\frac{1}{A}(a_2+b_2x+c_2y)v_2+\frac{1}{A}(a_3+b_3x+c_3y)v_3\\ &=N_1(x,y)v_1+N_2(x,y)v_2+N_3(x,y)v_3 \end{aligned}\tag{1-61}$
可以将二维三角形单元内任意一点的位移写成矩阵的形式，如下
$\begin{aligned} \mathbf u(x,y) = \begin{bmatrix} u(x,y)\\ v(x,y)\\ \end{bmatrix}&= \begin{bmatrix} N_1(x,y) & 0 & N_2(x,y) & 0 & N_3(x,y) & 0\\ 0 & N_1(x,y) & 0 & N_2(x,y) & 0 & N_3(x,y)\\ \end{bmatrix} \begin{bmatrix} u_1\\v_1\\u_2\\v_2\\u_3\\v_3 \end{bmatrix}\\&=\mathbf N(x,y)\cdot \mathbf q^e \end{aligned}\tag{1-62}$
我们称 $\mathbf N(x,y)$ 为形状函数矩阵，而 $\mathbf q^e$ 为单元的节点位移列阵。形状函数矩阵作用就是用节点位移列阵插值得到任意一点位移的插值函数。

1.2.2.3 单元应变

应变和位移的关系由几何方程确定，在本文二维的例子中，应变仅有三个分量，那么应变与位移的关系可以如下式所示。
$\begin{aligned} \boldsymbol{\varepsilon}(x,y) = \begin{bmatrix} \varepsilon_{xx}\\ \varepsilon_{yy}\\ \gamma_{xy} \end{bmatrix}= \begin{bmatrix} \frac{\partial u}{\partial x}\\ \frac{\partial v}{\partial y}\\ \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \end{bmatrix}&= \begin{bmatrix} \frac{\partial }{\partial x} & 0\\ 0 & \frac{\partial }{\partial y}\\ \frac{\partial }{\partial y} & \frac{\partial }{\partial x} \end{bmatrix}\cdot \begin{bmatrix} u(x,y)\\ v(x,y) \end{bmatrix}\\&=\begin{bmatrix} \frac{\partial }{\partial x} & 0\\ 0 & \frac{\partial }{\partial y}\\ \frac{\partial }{\partial y} & \frac{\partial }{\partial x} \end{bmatrix}\cdot\mathbf N(x,y)\cdot \mathbf q^e\\&=\mathbf B(x,y)\cdot \mathbf q^e \end{aligned}\tag{1-63}$

1.2.2.4 单元应力

应力应变关系又成本构方程，在线性范围内，应力与应变的关系写成矩阵如下(本例中为二维三角形单元，应力分量只有三个)
$\begin{aligned} \boldsymbol{\sigma}(x,y) = \begin{bmatrix} \sigma_{xx}\\ \sigma_{yy}\\ \tau_{xy} \end{bmatrix}&= \frac{E}{1-\mu^2}\begin{bmatrix} 1 & \mu & 0\\ \mu & 1 & 0\\ 0 & 0 & \frac{1-\mu}{2} \end{bmatrix}\begin{bmatrix} \varepsilon_{xx}\\ \varepsilon_{yy}\\ \gamma_{xy} \end{bmatrix}\\ &=\mathbf D(x,y)\cdot \boldsymbol{\varepsilon}(x,y)\\ &=\mathbf D(x,y)\cdot \mathbf B(x,y)\cdot \mathbf q^e\\ &=\mathbf S(x,y)\cdot \mathbf q^e \end{aligned}\tag{1-64}$

1.2.2.5 单元刚度矩阵

回顾物体的总势能表达式，并将其写成矩阵的形式
$\begin{aligned} II&=\frac{1}{2}\int_{\Omega}\sigma_{ij}\varepsilon_{ij}d\Omega-\int_{\Omega} b_{i}u_{i}d\Omega -\int_{A}p_{i}u_{i}dA\\ &=\frac{1}{2}\int_{\Omega}\boldsymbol{\sigma}^T\boldsymbol{\varepsilon}\mathbf d\Omega-\int_{\Omega} \mathbf b^T \mathbf u \mathbf d\Omega -\int_{A}\mathbf p^T\mathbf u \mathbf dA \end{aligned}\tag{1-65}$
将前述单元的应力、应变矩阵等代入上式，建立单元的总势能表达式，如下所示

$\begin{aligned} II &=\frac{1}{2}\int_{\Omega}\boldsymbol{\sigma}^T\boldsymbol{\varepsilon}\mathbf d\Omega-\int_{\Omega} \mathbf b^T \mathbf u \mathbf d\Omega -\int_{A}\mathbf p^T\mathbf u \mathbf dA\\ &=\frac{1}{2}\int_{\Omega}\boldsymbol{q^e}^T\boldsymbol{B^TD^TBq^e}\mathbf d\Omega-\int_{\Omega} \mathbf b^T \boldsymbol{ Nq^e d}\Omega -\int_{A}\mathbf p^T\boldsymbol{ Nq^e d}A\\ &=\frac{1}{2}\boldsymbol{q^e}^T(\int_{\Omega}\boldsymbol{B^TD^TB d}\Omega)\boldsymbol{q^e}-(\int_{\Omega} \mathbf b^T \boldsymbol{ Nd}\Omega +\int_{A}\mathbf p^T\boldsymbol{ Nd}A)\boldsymbol{q^e}\\ &=\frac{1}{2}\boldsymbol{q^e}^T\boldsymbol{K^e}\boldsymbol{q^e}-\boldsymbol{P^e}^T\boldsymbol{q^e} \end{aligned}\tag{1-66}$
回顾之前最小势能原理，真实位移是使得势能取得最小值，上式是单元的总势能，那么真实的位移使得势能取到最小值，也就是说总势能一阶变分为零，即
$\delta II = \frac{\partial II}{\partial q^e}\delta q^e=(\boldsymbol{K^e}\boldsymbol{q^e}-\boldsymbol{P^e})\delta q^e=0 \tag{1-67}$
由于 $\delta q^e$ 的任意性，可以得到
$\boldsymbol{K^e}\boldsymbol{q^e}-\boldsymbol{P^e}=0\tag{1-68}$
其中 $\boldsymbol{K^e}$ 为单元刚度矩阵（由上式中的量纲分析， $\boldsymbol{K^e}$ 为刚度的量纲）， $\boldsymbol{P^e}$ 为节点等效载荷。

1.2.2.6 整体刚度矩阵

对于物体来说，一般有多个单元，那么就可以建立多个单元平衡方程，如下所示
$\boldsymbol{K^e_{(i)}}\boldsymbol{q^e_{(i)}}-\boldsymbol{P^e_{(i)}}=0\tag{1-69}$
按照顺序将所有的单元的节点位移组成总体节点位移列阵，相应的节点等效载荷和单元刚度矩阵组成总体节点等效载荷列阵和整体刚度矩阵，以下图为例，下图总共有两个单元，得到两个单元平衡方程
$\begin{aligned} \begin{bmatrix} k^1_{11} & k^1_{12} & k^1_{13} & k^1_{14} & k^1_{15} & k^1_{16}\\ k^1_{21} & k^1_{22} & k^1_{23} & k^1_{24} & k^1_{25} & k^1_{26}\\ k^1_{31} & k^1_{32} & k^1_{33} & k^1_{34} & k^1_{35} & k^1_{36}\\ k^1_{41} & k^1_{42} & k^1_{43} & k^1_{44} & k^1_{45} & k^1_{46}\\ k^1_{51} & k^1_{52} & k^1_{53} & k^1_{54} & k^1_{55} & k^1_{56}\\ k^1_{61} & k^1_{62} & k^1_{63} & k^1_{64} & k^1_{65} & k^1_{66} \end{bmatrix}\begin{bmatrix} u_{1}\\ v_{1}\\u_{2}\\ v_{2}\\u_{3}\\ v_{3}\\ \end{bmatrix}=\begin{bmatrix} p^1_{x1}\\ p^1_{y1}\\p^1_{x2}\\ p^1_{y2}\\p^1_{x3}\\ p^1_{y3}\\ \end{bmatrix}\\ \begin{bmatrix} k^2_{11} & k^2_{12} & k^2_{13} & k^2_{14} & k^2_{15} & k^2_{16}\\ k^2_{21} & k^2_{22} & k^2_{23} & k^2_{24} & k^2_{25} & k^2_{26}\\ k^2_{31} & k^2_{32} & k^2_{33} & k^2_{34} & k^2_{35} & k^2_{36}\\ k^2_{41} & k^2_{42} & k^2_{43} & k^2_{44} & k^2_{45} & k^2_{46}\\ k^2_{51} & k^2_{52} & k^2_{53} & k^2_{54} & k^2_{55} & k^2_{56}\\ k^2_{61} & k^2_{62} & k^2_{63} & k^2_{64} & k^2_{65} & k^2_{66} \end{bmatrix}\begin{bmatrix} u_{2}\\ v_{2}\\u_{3}\\ v_{3}\\u_{4}\\ v_{4}\\ \end{bmatrix}=\begin{bmatrix} p^2_{x2}\\ p^2_{y2}\\p^2_{x3}\\ p^2_{y3}\\p^2_{x4}\\ p^2_{y4}\\ \end{bmatrix}\end{aligned}\tag{1-70}$
叠加整理后，整体平衡方程为
$\begin{bmatrix} k^1_{11} & k^1_{12} & k^1_{13} & k^1_{14} & k^1_{15} & k^1_{16} & 0 & 0\\ k^1_{21} & k^1_{22} & k^1_{23} & k^1_{24} & k^1_{25} & k^1_{26} & 0 & 0\\ k^1_{31} & k^1_{32} & k^1_{33}+k^2_{11} & k^1_{34}+k^2_{12} & k^1_{35}+k^2_{13} & k^1_{36}+k^2_{14} & k^2_{15} & k^2_{16}\\ k^1_{41}& k^1_{42} & k^1_{43}+k^2_{21}& k^1_{44}+k^2_{22}& k^1_{45} +k^2_{23} & k^1_{46}+k^2_{24} & k^2_{25} & k^2_{26}\\ k^1_{51}& k^1_{52}& k^1_{53}+k^2_{31} & k^1_{54}+k^2_{32} & k^1_{55}+k^2_{33} & k^1_{56}+k^2_{34} & k^2_{35} & k^2_{36}\\ k^1_{61}& k^1_{62}& k^1_{63}+k^2_{41} & k^1_{64}+k^2_{42} & k^1_{65}+k^2_{43} & k^1_{66}+k^2_{44} & k^2_{45} & k^2_{46}\\ 0 & 0 & k^2_{51} & k^2_{52} & k^2_{53} & k^2_{54} & k^2_{55} & k^2_{56}\\ 0 & 0 & k^2_{61} & k^2_{62} & k^2_{63} & k^2_{64} & k^2_{65} & k^2_{66} \end{bmatrix}\begin{bmatrix} u_{1}\\ v_{1}\\u_{2}\\ v_{2}\\u_{3}\\ v_{3}\\u_{4}\\ v_{4}\\ \end{bmatrix}=\begin{bmatrix} p^1_{x1}\\ p^1_{y1}\\p^1_{x2}+p^2_{x2}\\ p^1_{y2}+p^2_{y2}\\p^1_{x3}+p^2_{x3}\\ p^1_{y3}+p^2_{y3}\\p^2_{x4}\\ p^2_{y4}\\ \end{bmatrix}\tag{1-71}$
在这里插入图片描述

1.2.2.7 处理约束

如果在上例中，假设节点1为约束，那么 $u_1=v_1=0$ ，上式变为
$\begin{bmatrix} k^1_{11} & k^1_{12} & k^1_{13} & k^1_{14} & k^1_{15} & k^1_{16} & 0 & 0\\ k^1_{21} & k^1_{22} & k^1_{23} & k^1_{24} & k^1_{25} & k^1_{26} & 0 & 0\\ k^1_{31} & k^1_{32} & k^1_{33}+k^2_{11} & k^1_{34}+k^2_{12} & k^1_{35}+k^2_{13} & k^1_{36}+k^2_{14} & k^2_{15} & k^2_{16}\\ k^1_{41}& k^1_{42} & k^1_{43}+k^2_{21}& k^1_{44}+k^2_{22}& k^1_{45} +k^2_{23} & k^1_{46}+k^2_{24} & k^2_{25} & k^2_{26}\\ k^1_{51}& k^1_{52}& k^1_{53}+k^2_{31} & k^1_{54}+k^2_{32} & k^1_{55}+k^2_{33} & k^1_{56}+k^2_{34} & k^2_{35} & k^2_{36}\\ k^1_{61}& k^1_{62}& k^1_{63}+k^2_{41} & k^1_{64}+k^2_{42} & k^1_{65}+k^2_{43} & k^1_{66}+k^2_{44} & k^2_{45} & k^2_{46}\\ 0 & 0 & k^2_{51} & k^2_{52} & k^2_{53} & k^2_{54} & k^2_{55} & k^2_{56}\\ 0 & 0 & k^2_{61} & k^2_{62} & k^2_{63} & k^2_{64} & k^2_{65} & k^2_{66} \end{bmatrix}\begin{bmatrix} 0\\ 0\\u_{2}\\ v_{2}\\u_{3}\\ v_{3}\\u_{4}\\ v_{4}\\ \end{bmatrix}=\begin{bmatrix} R_{x1}\\ R_{y1}\\p^1_{x2}+p^2_{x2}\\ p^1_{y2}+p^2_{y2}\\p^1_{x3}+p^2_{x3}\\ p^1_{y3}+p^2_{y3}\\p^2_{x4}\\ p^2_{y4}\\ \end{bmatrix}\tag{1-72}$
其中 $R_{x1}、R_{y1}$ 为支反力，将上式化为两个方程
$\begin{aligned}\begin{bmatrix} k^1_{33}+k^2_{11} & k^1_{34}+k^2_{12} & k^1_{35}+k^2_{13} & k^1_{36}+k^2_{14} & k^2_{15} & k^2_{16}\\ k^1_{43}+k^2_{21}& k^1_{44}+k^2_{22}& k^1_{45} +k^2_{23} & k^1_{46}+k^2_{24} & k^2_{25} & k^2_{26}\\ k^1_{53}+k^2_{31} & k^1_{54}+k^2_{32} & k^1_{55}+k^2_{33} & k^1_{56}+k^2_{34} & k^2_{35} & k^2_{36}\\ k^1_{63}+k^2_{41} & k^1_{64}+k^2_{42} & k^1_{65}+k^2_{43} & k^1_{66}+k^2_{44} & k^2_{45} & k^2_{46}\\ k^2_{51} & k^2_{52} & k^2_{53} & k^2_{54} & k^2_{55} & k^2_{56}\\ k^2_{61} & k^2_{62} & k^2_{63} & k^2_{64} & k^2_{65} & k^2_{66} \end{bmatrix}\begin{bmatrix} u_{2}\\ v_{2}\\u_{3}\\ v_{3}\\u_{4}\\ v_{4}\\ \end{bmatrix}&=\begin{bmatrix} p^1_{x2}+p^2_{x2}\\ p^1_{y2}+p^2_{y2}\\p^1_{x3}+p^2_{x3}\\ p^1_{y3}+p^2_{y3}\\p^2_{x4}\\ p^2_{y4}\\ \end{bmatrix}\\ \overline K \overline q^e&=\overline P \end{aligned}\tag{1-73}$

$\begin{bmatrix} k^1_{13} & k^1_{14} & k^1_{15} & k^1_{16}\\ k^1_{23} & k^1_{24} & k^1_{25} & k^1_{26}\\ \end{bmatrix}\begin{bmatrix} u_{2}\\ v_{2}\\u_{3}\\ v_{3} \end{bmatrix}=\begin{bmatrix} R_{x1}\\ R_{y1}\\ \end{bmatrix}\tag{1-74}$

1.2.2.8 求解节点载荷列阵

1.2.2.9 求解位移列阵

最终通过上式求得位移列阵 $\overline q^e=\overline K^{-1}\cdot \overline P$ ，再将与 $u_1,v_1]^T=[0,0]^T$ 组装成最终的节点位移列阵 $\boldsymbol{q^e}$ 。

1.2.2.10 计算应力矩阵等

计算得到最终的节点位移列阵 $\boldsymbol{q^e}$ 后通过一系列的回代得到单元内任意位置的位移、应变、应力矩阵，而式（1-74）用来求解约束支反力。
$\begin{aligned} \mathbf u(x,y) =\mathbf N(x,y)\cdot \mathbf q^e\\ \boldsymbol{\varepsilon}(x,y) = \mathbf B(x,y)\cdot \mathbf q^e\\ \boldsymbol{\sigma}(x,y) = \mathbf S(x,y)\cdot \mathbf q^e \end{aligned}\tag{1-75}$