🌟有限元分析学习笔记（一）：刚度矩阵的推导与组装

完成时间：2021.12.28
感谢 江南大学胡清元 老师的讲解
参考：

• 《有限单元法》-- 王勖成

• 平面四边形单元有限元实现——有限元实践笔记（6）

• 有限元分析：ANSYS理论与应用 第四版

• 四节点矩形单元有限元分析过程

前言

本人是信息与计算科学专业的大四学生，因工作和毕业论文需要，所以学习有限元，本人除了高中物理以外没有任何其他力学基础，所以这是一篇纯数学和程序员向的有限元笔记，与网上其他教程不同，除去了有限元中繁冗复杂的力学方程和物理条件，只剩下数学推导和少部分的物理方程，不妥之处可以发送内容到邮箱：hi@qiql.net, 阅读本文需要你有以下学科知识：微积分、线代、数值分析 或了解以下知识点：求导、微分、链式法则、矩阵乘法、求逆、插值、二维积分、雅可比矩阵、数值积分(高斯积分)。

假如读者与我一样，无力学背景，建议读者在掌握了上述前两个知识点的情况下，先看 有限元分析：ANSYS理论与应用 第四版 这本书中的绪论部分中的第一个例子（必须能完全理解该例子，物理过程不难，有初中物理基础即可），了解一下有限元的基本思想；然后阅读《有限单元法》-- 王勖成 书中的第一章和第二章，这本书若无力学基础，读来十分晦涩，可连滚带爬往前阅读，无需过多在意其中复杂的物理方程；然后带着无尽的疑惑阅读四节点矩形单元有限元分析过程，该文为PPT，可能读者访问时链接已失效，可以向我邮件求助，该PPT讲的十分清晰，但其中的雅克比变换有些不符合常规；理清思路后，结合本文和平面四边形单元有限元实现——有限元实践笔记（6） 一起阅读，慢慢理解，注意：本文与后者中的符号有些许不同，十分容易造成阅读混乱，逻辑不清晰，读者需要仔细认真推导其中方程，不断搜索与查找，不能只看不写，哪怕抄一抄也是极好的，且后者中的雅克比矩阵部分有一些错误，推导时注意规避。

一、设位移函数

考虑如图 1 的Cook短斜梁问题：

图1 Cook 短斜梁
对其进行离散化，假如划分后的网格中的某单元如图 2 所示：

假如单元节点序号如图2所示：

图2 Cook 短斜梁中的一个单元

对于平面上的四节点单元，每个结点有两个自由度

$$\bold a_i= \left( \begin{array}{ccc} \bold u_i\\ \bold v_i\\ \end{array}\tag{1.0} \right)$$

即一共有 8 个结点位移

$$\bold a^e= \left( \begin{array}{ccc} \bold a_i\\ \bold a_j\\ \bold a_m\\ \bold a_n\\ \end{array} \right) =\begin{bmatrix} u_i & v_i & u_j & v_j & u_m & v_m & u_n & v_n\\ \end{bmatrix}^T$$

设位移函数为：

$$\begin{cases} u = \beta_1+\beta_2\cdot x + \beta_3\cdot y+\beta_4\cdot x \cdot y\\ v = \beta_5+\beta_6\cdot x + \beta_7\cdot y+\beta_8\cdot x \cdot y\\ \end{cases}$$

或者，写为：

$$\{\delta(x,y)\}= \begin{bmatrix} f(x,y)\\ \end{bmatrix} \cdot \{\bold\beta\}$$

整理，得到：

$$\begin{bmatrix} u_i\\ v_i\\ u_j\\ v_j\\ u_m\\ v_m\\ u_n\\ v_n\\ \end{bmatrix} = \underbrace{ \begin{bmatrix} 1 &x_i & y_i & x_i\cdot y_i & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 &x_i & y_i & x_i\cdot y_i\\ 1 &x_j & y_j & x_j\cdot y_j & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 &x_j & y_j & x_j\cdot y_j\\ 1 &x_m & y_m & x_m\cdot y_m & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 &x_m & y_m & x_m\cdot y_m\\ 1 &x_n & y_n & x_n\cdot y_n & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 &x_n & y_n & x_n\cdot y_n\\ \end{bmatrix} }_{设为\bold A_{(8*8)}} \cdot \begin{bmatrix} \beta_1\\ \beta_2\\ \beta_3\\ \beta_4\\ \beta_5\\ \beta_6\\ \beta_7\\ \beta_8 \end{bmatrix}$$

二、求解位移函数中的未知量

得到了 $\bold a^e=\bold A \cdot \beta$ 后解之，得

$$\{\beta\}= \bold A ^{-1} \cdot \bold a^e \tag 6$$

将（6）式代入（4）中，得到：

$$\{\delta(x,y)\}= \begin{bmatrix} f(x,y)\\ \end{bmatrix} \cdot \{\bold\beta\} = \begin{bmatrix} f(x,y)\\ \end{bmatrix} \underbrace{ \cdot \bold A^{-1}\cdot \bold a^e }_{ \{\beta\}} =[N(x,y)]\cdot \bold a^e \tag 7$$

整理，得到：

$$\{\delta(x,y)\}= \begin{bmatrix} u\\ v \end{bmatrix} = \begin{bmatrix} N_i & 0 & N_j & 0 & N_m & 0 & N_n & 0\\ 0 & N_i & 0 & N_j & 0 & N_m & 0 & N_n \end{bmatrix} \cdot \begin{bmatrix} u_i\\ v_i\\ u_j\\ v_j\\ u_m\\ v_m\\ u_n\\ v_n\\ \end{bmatrix}$$

其中，N(x,y) 就是 形函数，或称之为 插值函数

三、形函数的性质

• 性质一

$$N(x_i,y_i)=\delta_{ij}= \begin{cases} 1&,j=i\\ 0&,j\neq i \end{cases}$$

• 性质二

$$\sum_i N_i(x,y)=N_i(x,y)+N_j(x,y)+N_m(x,y)+N_n(x,y) = 1$$

​ 即，在单元内的任意一个位置，形函数的和均为 1

在本文中，因为单元是四个结点

假如，四个结点的坐标为：

结点序号	x坐标	y坐标
$i$	$-a$	$-b$
$j$	$a$	$-b$
$m$	$c$	$d$
$n$	$-c$	$d$

表1，四节点单元的坐标

那么，形函数为：

$$\left\{ \begin{aligned} N_i(x,y) & = \frac{1}{4} \cdot (1-\frac{x}{a})\cdot (1-\frac{y}{b})\\ N_j(x,y) &= \frac{1}{4} \cdot (1+\frac{x}{a})\cdot (1-\frac{y}{b})\\ N_m(x,y) &= \frac{1}{4} \cdot (1+\frac{x}{c})\cdot (1+\frac{y}{d})\\ N_n(x,y) &= \frac{1}{4} \cdot (1-\frac{x}{c})\cdot (1+\frac{y}{d}) \end{aligned} \right.$$

至此，整理，得到：

$$\begin{aligned} & u=N_i(x,y)\cdot u_i + N_j(x,y) \cdot u_j + N_m(x,y) \cdot u_m + N_n(x,y) \cdot u_n\\ & v=N_i(x,y)\cdot v_i + N_j(x,y) \cdot v_j + N_m(x,y) \cdot v_m + N_n(x,y) \cdot v_n \end{aligned}$$

因为式中的 $u_i,v_i(i,j,m,n)$ 是常数，所以位移函数是完全由形函数来决定的

对得到的形函数，验证形函数的性质（1）和性质（2）

• 验证性质（1）： $N(x_i,y_i)=\delta_{ij}= \begin{cases} 1&,j=i\\ 0&,j\neq i \end{cases}$

$$\left\{ \begin{aligned} N_i(x_i,y_i)&=\frac{1}{4} \cdot (1-\frac{-a}{a})\cdot (1-\frac{-b}{b}) = 1\\ N_i(x_j,y_j) &=\frac{1}{4} \cdot (1-\frac{a}{a})\cdot (1-\frac{-b}{b}) = 0\\ N_i(x_m,y_m) &=\frac{1}{4} \cdot (1-\frac{-c}{a})\cdot (1-\frac{d}{d}) = 0\\ N_i(x_n,y_n) &=\frac{1}{4} \cdot (1-\frac{c}{a})\cdot (1-\frac{d}{d}) = 0 \end{aligned} \right.$$

对于结点 $j,m,n$ 同理

• 验证性质（2）： $\sum_i N_i(x,y)=N_i+N_j(x,y)+N_m(x,y)+N_n(x,y) = 1$

$$\begin{aligned} \sum_i N_i(x,y) & =N_i+N_j(x,y)+N_m(x,y)+N_n(x,y)\\ &=\frac{1}{4} \cdot (1-\frac{x}{a})\cdot (1-\frac{y}{b})+ \frac{1}{4} \cdot (1+\frac{x}{a})\cdot (1-\frac{y}{b})+ \frac{1}{4} \cdot (1+\frac{x}{c})\cdot (1+\frac{y}{b})+ \frac{1}{4} \cdot (1-\frac{x}{c})\cdot (1+\frac{y}{b})\\ & = \frac{1}{4} \cdot[ (1-\frac{y}{b})(1-\frac{x}{a}+1+\frac{x}{a})+(1-\frac{y}{b})(1+\frac{x}{c}+1-\frac{x}{c})]\\ & = \frac{1}{4} \cdot [(1-\frac{y}{b})\cdot2 + (1+\frac{y}{b}\cdot 2)\\ & =1 \end{aligned}$$

为方便后文叙述，以及符号规范，令 $\bold u = \begin{bmatrix}u\\ v \end{bmatrix}$ 即：用 $\bold u$ 表示位移

综上，得到：

$$\bold u = \bold N \cdot \bold a^e$$

其中：

$$\bold N = \begin{bmatrix} N_i(x,y) & 0 & N_j(x,y) & 0 & N_m(x,y) & 0 & N_n(x,y) & 0\\ 0 & N_i(x,y) & 0 & N_j(x,y) & 0 & N_m(x,y) & 0 & N_n(x,y)\\ \end{bmatrix}$$

四、根据几何方程求应变矩阵

几何方程的矩阵形式为：

$$\bold \epsilon =\bold L\bold u$$

其中， $\bold L$ 为微分算子， $\bold u$ 为位移，在二维条件下，其值为：

$$\bold L =\begin{bmatrix} \frac{\partial}{\partial x} & 0\\ 0 & \frac{\partial}{\partial y}\\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y}\\ \end{bmatrix} =A^T$$

若是在三维条件下，其值为：

$$\bold L =\begin{bmatrix} \frac{\partial}{\partial x} & 0 & 0\\ 0 & \frac{\partial}{\partial y} & 0\\ 0 & 0 & \frac{\partial}{\partial z}\\ \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0\\ 0 & \frac{\partial}{\partial z} & \frac{\partial}{\partial y}\\ \frac{\partial}{\partial z} & 0 & \frac{\partial}{\partial x}\\ \end{bmatrix} =A^T$$

根据（11）式、（13）式、（14）式，本文所研究问题（二维平面应力问题）的几何方程为：

$$\{\bold \epsilon \}= \begin{pmatrix} \epsilon_x \\ \epsilon_y \\ \gamma_{xy} \end{pmatrix} = \underbrace{ \begin{bmatrix} \frac{\partial}{\partial x} & 0\\ 0 & \frac{\partial}{\partial y}\\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y}\\ \end{bmatrix} \cdot \bold N }_{\bold B} \cdot \bold a^e =\bold B\cdot \bold a^e$$

其中：

• $\gamma_{xy}$ 是剪应变 $\epsilon_x$ 和 $\epsilon_y$ 分别是 $x$ 和 $y$ 方向的正应变
• $\bold N$ 为形函数矩阵（2*8）, 其值为：

$$\bold N = \begin{bmatrix} N_i & 0 & N_j & 0 & N_m & 0 & N_n & 0\\ 0 & N_i & 0 & N_j & 0 & N_m & 0 & N_n \end{bmatrix}$$

将其进行分块，不妨写为：

$$\begin{aligned} & \begin{bmatrix} \bold N_i & \bold N_j &\bold N_m &\bold N_n\\ \end{bmatrix} & 其中，\bold N_i= \begin{bmatrix} N_i & 0\\ 0 & N_i \end{bmatrix} & ,\;\;\; (i,j,m,n) \end{aligned}$$

其中，

• $\bold a^e$ 为单元位移列阵， $\bold B$ 即为应变矩阵，对形函数矩阵 $\bold N$ 求微分，即可得到 $\bold B$
• 显然，对 $\bold N$ 进行分块后，对 $\bold B$ 也可以进行分块，不妨写为：

$$\begin{aligned} & \bold B = \bold L \cdot \bold N = \begin{bmatrix} \bold B_i & \bold B_j &\bold B_m &\bold B_n\\ \end{bmatrix} & 其中，\bold B_i= \begin{bmatrix} \frac{\partial}{\partial x} & 0\\ 0 & \frac{\partial}{\partial y}\\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y}\\ \end{bmatrix} \begin{bmatrix} N_i & 0\\ 0 & N_i \end{bmatrix} =\bold L \cdot \bold N_i & \;,\;\;(i,j,m,n) \end{aligned}$$

五、根据物理方程求应力矩阵

弹性力学中，将应力与应变之间的转换关系也称为弹性关系。对于各向同性的线弹性材料，应力通过应变的表达式可用矩阵形式表示为：

$$\bold \sigma = \bold D \cdot \epsilon$$

其中，$\epsilon$ 为应变矩阵，$\sigma$ 为应力矩阵，D 为弹性矩阵

$$\bold D =\frac{E \cdot (1-\nu)}{(1+\nu)\cdot (1-2\nu)} =\begin{bmatrix} 1 & \frac{\nu}{(1-\nu)} & \frac{\nu}{(1-\nu)} & 0 & 0 & 0\\ \frac{\nu}{(1-\nu)} & 1 & \frac{\nu}{(1-\nu)} & 0 & 0 & 0\\ \frac{\nu}{(1-\nu)} & \frac{\nu}{(1-\nu)} & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} & 0 & 0\\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} & 0\\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} \end{bmatrix}$$

如你所见，$\bold D$ 是一个对称阵，作为弹性矩阵，它完全取决于弹性体材料的弹性模量 $E$ 和泊桑比 $\nu$ 。

另外，此矩阵还有另一种表示方式，也可以采用拉梅常数 $G$ 和 $\lambda$ ，它们和 $E$，$\nu$ 的关系如下

$$\begin{aligned} &G=\frac{E}{2(1+\nu)}\;\;\; &\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)} \end{aligned}$$

$G$ 也称为剪切弹性模量，注意到：

$$\lambda + 2 G = \frac {E(1-\nu)}{(1+\nu)(1-2\nu)}$$

故弹性矩阵 $\bold D$ 亦可以表示为：

$$\bold D =\begin{bmatrix} \lambda + 2G & \lambda & \lambda & 0 & 0 & 0\\ \lambda & \lambda + 2G & \lambda & 0 & 0 & 0\\ \lambda &\lambda & \lambda + 2G & 0 & 0 & 0\\ 0 & 0 & 0 & G & 0 & 0\\ 0 & 0 & 0 & 0 & G & 0\\ 0 & 0 & 0 & 0 & 0 & G\\ \end{bmatrix}$$

所以，物理方程的另一种形式是：

$$\epsilon = \bold C \sigma$$

其中， $\bold C$ 是柔度矩阵，它和弹性矩阵是互逆关系

注意：本文中在求解应力矩阵时，暂不采用柔度矩阵的方式

根据《有限单元法》一书中第 61 页的说明：

• 对于平面应力问题：

$$\begin{aligned} & E_0=E & \nu_0 = \nu \end{aligned}$$

• 对于平面应变问题：

$$\begin{aligned} & E_0=\frac{E}{1-\nu^2} & \nu_0=\frac{v}{1-\nu} \end{aligned}$$

综上，求得 $\bold D$ 矩阵为：

$$\bold D = \frac{E_0}{1-\nu_0^2}\cdot \begin{bmatrix} 1 & \nu_0 & 0 \\ \nu_0 & 1 &0 \\ 0 & 0 & \frac{1-\nu_0}{2} \end{bmatrix}$$

那么：

$$\{\sigma\} = [\bold D]\cdot\{\epsilon\} = \underbrace{ [\bold D]\cdot [\bold B] }_{[\bold S]} \cdot\{\bold a^e\}=[\bold S]\cdot \{\bold a^e\}$$

其中， $\bold S$ 即为应力矩阵

由于， $\bold B$ 可以写为分块矩阵，所以 $\bold S$ 也可以写为下面形式的分块矩阵：

$$\begin{aligned} & \bold S = \bold D \cdot \bold B= \begin{bmatrix} \bold S_i & \bold S_j &\bold S_m &\bold S_n\\ \end{bmatrix} & 其中，\bold S_i =\bold D \cdot \bold B_i & \;,\;\;(i,j,m,n) \end{aligned}$$

注意： $\bold S$ 与 $\bold B$ 相同，也都是常量阵，但在很多情况下，不单独定义应力矩阵 $\bold S$ ，而直接用 $\bold D\bold B$ 去计算应力

六、转换单元坐标系

虚功原理这块比较复杂，但结论是：

$$\bold K^e=\iint\limits_{A}\bold B^T\cdot \bold D \cdot \bold B \cdot t \;\;\mathrm{d}x\mathrm{d}y$$

其中， $t$ 为单元的厚度，$\bold K^e$ 是单元刚度阵

如果要直接在单元内进行积分，由于各个单元的形状、大小不仅相同，求解起来变得复杂，而且工程问题也不一定非要取得精确解，有误差可控的数值解就可以了，如果能将积分区域规定到一个固定的区域内，就可以采用很多现有的数值积分方法来进行快速求解，具体的说，就是致力于寻找一个映射，将原来的单元，无论形状如何，都可以映射到一个基准单元中，积分区域将由原来的单元转变为基准单元。（这段是我自己的理解）

在正式推导之前，我们需要更改原来的形函数，之前在划分单元后，每个单元的形状都不尽相同，所以位移插值函数也不同，求解效率低，通过将形状不规整的实际单元和形状规整的标准单元之间建立坐标映射关系，在标准单元上去研究它的性质，却不改变原来的问题，映射称为参数映射，单元称为参数单元，如果坐标变换和试函数里的形函数、插值节点完全相同，则变换就称为等参数变换，参数单元称等参数单元。相应的还有亚参元和超参元。^[1]

一种粗暴的变换方式是直接进行压缩，如这篇：https://www.docin.com/p-764314942.html

但是我感觉这样映射其实没啥用，除了去掉了量纲，一种更常见的映射如图 3 所示：

图3 自然坐标系到基准坐标系下的映射示例图

这一点在《有限单元法》-- 王勖成 第二章中尚未提到，是胡老师讲过之后我才知道，顺带搜了一下，其思想在于：让基准单元到每个物理单元之间存在一个映射，比较容易想到的是就是线性映射了，在那篇知乎文章中，映射是这样的：

$$\begin{aligned} & x=\frac{1}{4}[(1-\xi)(1-\eta)x_i+(1+\xi)(1-\eta)x_j+(1+\xi)(1+\eta)x_m+(1-\xi)(1+\eta)x_n]\\ & y=\frac{1}{4}[(1-\xi)(1-\eta)y_i+(1+\xi)(1-\eta)y_j+(1+\xi)(1+\eta)y_m+(1-\xi)(1+\eta)y_n] \end{aligned}$$

根据这个映射，可以将母单元上的四个点一一映射到物理单元上的四个点，注意：$x_i \;,\;(i,j,m,n)$ 是已知的，所以根据母单元的中任意一点坐标 ($\xi,\eta$) 可以求得其在物理单元中的坐标 ($x_0,y_0$) ，反过来也一样。

更进一步，将公式（32）写为矩阵形式：

$$\begin{bmatrix} x\\ y \end{bmatrix} = \begin{bmatrix} N_i & 0 & N_j & 0 & N_m & 0 & N_n & 0\\ 0 & N_i & 0 & N_j & 0 & N_m & 0 & N_n \end{bmatrix} \cdot \begin{bmatrix} x_i\\ y_i\\ x_j\\ y_j\\ x_m\\ y_m\\ x_n\\ y_n\\ \end{bmatrix}$$

容易得到，形函数为：

$$\left\{ \begin{aligned} N_i(\xi,\;\eta) &=\frac{1}{4} \cdot (1-\xi)\cdot (1-\eta)\\ N_j(\xi,\;\eta) &=\frac{1}{4} \cdot (1+\xi)\cdot (1-\eta)\\ N_m(\xi,\;\eta) &=\frac{1}{4} \cdot (1+\xi)\cdot (1+\eta)\\ N_n(\xi,\;\eta) &=\frac{1}{4} \cdot (1-\xi)\cdot (1+\eta)\\ \end{aligned} \right.$$

那么，因为之前推导的方程是在原坐标系下，所以现在需要转换到基准坐标系下来

我们当然可以直接求 一个变换 $f$ ，使得：

$$N_i(\xi,\eta) = f\cdot N_i(x,y)$$

但是考虑到：

$$\begin{aligned} & \bold B = \bold L \cdot \bold N & \bold S = \bold D \cdot B = \bold D \cdot L \cdot \bold N \end{aligned}$$

所以不妨更近一步，考虑新 $\bold B$ 与旧 $\bold B$ 矩阵之间的转换

回到方程（16）

$$\{\bold \epsilon \}= \begin{pmatrix} \epsilon_x \\ \epsilon_y \\ \gamma_{xy} \end{pmatrix}= \bold L \cdot \bold u= \underbrace{ \begin{bmatrix} \frac{\partial}{\partial x} & 0\\ 0 & \frac{\partial}{\partial y}\\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y}\\ \end{bmatrix} \cdot \bold N }_{\bold B} \cdot \bold a^e =\bold B\cdot \bold a^e$$

由于坐标系发生了变化，所以我们现在需要求在基准坐标系下的应变

根据链式法则：

$$\frac{\partial f}{\partial \xi} = \frac{\partial f}{\partial x} \cdot \frac{\partial x}{\partial \xi} + \frac{\partial f}{\partial y} \cdot \frac{\partial y}{\partial \xi}\\ \frac{\partial f}{\partial \eta} = \frac{\partial f}{\partial x} \cdot \frac{\partial x}{\partial \eta} + \frac{\partial f}{\partial y} \cdot \frac{\partial y}{\partial \eta}\\$$

求解上式得到：

$$\begin{aligned} & \frac{\partial}{\partial x} = \frac{1}{\det[\bold J]}\left[ \frac{\partial x}{\partial \eta}\cdot\frac{\partial}{\partial \xi} - \frac{\partial x}{\partial \xi}\cdot\frac{\partial}{\partial \eta} \right]\\ & \frac{\partial}{\partial y} = \frac{1}{\det[\bold J]}\left[ \frac{\partial y}{\partial \xi}\cdot\frac{\partial}{\partial \eta} - \frac{\partial y}{\partial \eta}\cdot\frac{\partial}{\partial \xi} \right]\\ \end{aligned}$$

其中 $\bold J$ 为雅克比矩阵，关于对雅克比矩阵的具体描述，请参见：https://www.zhihu.com/question/22739036/answer/442131829

在本文中：

$$\begin{cases} \bold J = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi}\\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{bmatrix} \\ \det [J] =\frac{1}{8}\bold X_c\cdot \begin{bmatrix} 0 & 1-\eta & \eta - \xi & \xi - 1\\ \eta-1 & 0 & \xi+1 & -\xi-1\\ \xi-\eta & -\xi - 1 & 0 & \eta +1\\ 1-\xi & \xi + 1 & -\eta - 1 & 0 \end{bmatrix} \cdot Y_c\;, & 其中:\bold X_c = \begin{pmatrix} x_i\\ x_j\\ x_m\\ x_n \end{pmatrix}^T , \bold Y_c = \begin{pmatrix} y_i\\ y_j\\ y_m\\ Y_n \end{pmatrix} \end{cases}$$

如此，便将原本对 $(x, y)$ 坐标系下的微分算子，通过雅克比矩阵转换到 $(ξ, η)$ 坐标系下了

新的微分算子为：

$$\bold L^{'} = \frac{1}{\det[\bold J]}\cdot \begin{bmatrix} \frac{\partial x}{\partial \eta}\cdot\frac{\partial}{\partial \xi} - \frac{\partial x}{\partial \xi}\cdot\frac{\partial}{\partial \eta} & 0 \\ 0 & \frac{\partial y}{\partial \xi}\cdot\frac{\partial}{\partial \eta} - \frac{\partial y}{\partial \eta}\cdot\frac{\partial}{\partial \xi}\\ \frac{\partial x}{\partial \eta}\cdot\frac{\partial}{\partial \xi} - \frac{\partial x}{\partial \xi}\cdot\frac{\partial}{\partial \eta} & \frac{\partial y}{\partial \xi}\cdot\frac{\partial}{\partial \eta} - \frac{\partial y}{\partial \eta}\cdot\frac{\partial}{\partial \xi} \end{bmatrix}$$

那么，在新的坐标系下：

$$\epsilon = \bold L^{'} \cdot \bold u\\ \bold B(\xi,\eta)= \bold L^{'} \cdot \bold N(\xi,\eta)= \frac{1}{\det[\bold J]}\cdot \begin{bmatrix} \bold B_i & \bold B_j & \bold B_m & \bold B_n \end{bmatrix}$$

七、求解单元刚度矩阵

因为已经转换了坐标系，所以之前的求解单元刚度矩阵的公式就可以更新为：

$$\bold K^e=\iint\limits_{A}\bold B(x,y)^T\cdot \bold D \cdot \bold B(x,y) \cdot t \;\;\mathrm{d}x\mathrm{d}y =\int_{-1}^{1}\int_{-1}^{1}\bold B(\xi,\eta)^T\cdot \bold D \cdot \bold B(\xi,\eta) \cdot t\cdot \left|\bold J \right| \;\mathrm{d}\xi\mathrm{d}\eta$$

同上， $t$ 为单元厚度

要求解此积分，精确解难以获得，可以通过数值积分获取，比较常用的一种数值积分方式为高斯积分

对于如下函数的积分：

$$I = \int_{-1}^{1}y\;\mathrm{d}x$$

假如 $y_i$ 为样本点， $W_i$ 为样本点对应的权值，$n$ 为样本点的个数，可将上式积分写为：

$$I = \int_{-1}^{1}y\;\mathrm{d}x =\sum_{i=1}^{n}W_i\cdot y_i$$

对于二维高斯积分，同理：

$$\begin{aligned} I &= \int_{-1}^{1}\int_{-1}^{1}f(\xi,\eta)\;\mathrm{d}\xi\mathrm{d}\eta\\ &=\int_{-1}^{1}\sum_{i}W_i\cdot f(\xi_i,\eta)\mathrm{d}\eta\\ &=\sum_{j}W_j\left[\sum_{i}W_i\cdot f(\xi_i,\eta_j)\right]\\ &=\sum_{j}\sum_{i}W_iW_j\cdot f(\xi_i,\eta_j)\\ &\end{aligned}$$

查表，可得部分高斯积分样本值和权值

样本点数	样本点值	权值
$1$	$x_1 = 0$	$2$
$2$	$x_1,x_2=±0.577 350 269 189 626$	$1$
$3$	$x_1,x_3=±0.774 596 669 241 483$ $x_2=0$	$5/9$ $8/9$
$4$	$x_1,x_4=±0.861 136 311 594 053$ $x_2,x_3=±0.339 981 043 584 856$	$0.347 854 845 137 454$ $0.652 145 154 862 546$

由于本文是二维积分问题，选四点即可。已知 $x_1 = 0.577 350 269 189 626$ , $x_2 = -0.577 350 269 189 626$

令：

结点序号	坐标值	坐标值	权值
$1$	$ξ_1 = x_2$	$η_1 = x_2$	$W_1=1*1=1$
$2$	$ξ_2 = x_1$	$η_2 = x_2$	$W_2=1*1=1$
$3$	$ξ_3 = x_1$	$η_3 = x_1$	$W_3=1*1=1$
$4$	$ξ_4 = x_2$	$η_4 = x_2$	$W_4=1*1=1$

根据高斯积分，便可将积分问题转变为求和问题：

$$\bold K^e=t\cdot \left|\bold J \right|\cdot\{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\\ \;\;\;\bold B(\xi_1,\eta_1)^T\cdot \bold D \bold \cdot B(\xi_1,\eta_1)\cdot W_1\\ +\bold B(\xi_2,\eta_2)^T\cdot \bold D \bold \cdot B(\xi_2,\eta_2)\cdot W_2\\ +\bold B(\xi_3,\eta_3)^T\cdot \bold D \bold \cdot B(\xi_3,\eta_3)\cdot W_3\\ +\bold B(\xi_4,\eta_4)^T\cdot \bold D \bold \cdot B(\xi_4,\eta_4)\cdot W_4\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\}$$

八、组装总刚度矩阵

在进行组装之前，我们需要先找到一个变换 $\bold G$ ，能够将单元刚度阵转换到总刚矩阵中。假如该结构的总结点数为 $l$

不难想到，$\bold G$ 的作用是提高阶数，在《有限单元法》中，对其的定义为：

$$\bold a^e = \bold G \bold a^t$$

其中， $\bold a^e$ 为结点自由度（8*1），如果忘了可以翻看本文的最开始那部分， $\bold a^t$ 是整个结构的结点位移列阵（ $2l * 1$），通过此方程，便可以使得单元的结点自由度数和结构的结点自由度数相同。
其中：

$$\bold G\;\; = \begin{bmatrix} 0&0&\cdots & 1 & 0 & \cdots & 0 & 0 & \cdots & 0 & 0& \cdots & 0 & 0& \cdots & 0\\ 0&0&\cdots & 0 & 1 & \cdots & 0 & 0 & \cdots & 0 & 0& \cdots & 0 & 0& \cdots & 0\\ 0&0&\cdots & 0 & 0 & \cdots & 1 & 0 & \cdots & 0 & 0& \cdots & 0 & 0& \cdots & 0\\ 0&0&\cdots & 0 & 0 & \cdots & 0 & 1 & \cdots & 0 & 0& \cdots & 0 & 0& \cdots & 0\\ 0&0&\cdots & 0 & 0 & \cdots & 0 & 0 & \cdots & 1 & 0& \cdots & 0 & 0& \cdots & 0\\ 0&0&\cdots & 0 & 0 & \cdots & 0 & 0 & \cdots & 0 & 1& \cdots & 0 & 0& \cdots & 0\\ 0&0&\cdots & 0 & 0 & \cdots & 0 & 0 & \cdots & 0 & 0& \cdots & 1 & 0& \cdots & 0\\ 0&0&\cdots & 0 & 0 & \cdots & 0 & 0 & \cdots & 0 & 0& \cdots & 0 & 1& \cdots & 0\\ \end{bmatrix}$$

同样的，可以得到：

$$\bold a^t = \bold G^T \bold a^e$$

将上式展开写一下：

$$\bold a^t = \bold G^T \bold a^e = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{bmatrix} \begin{bmatrix} u_i \\ v_i \\ u_j \\ v_j \\ u_m \\ v_m \\ u_n \\ v_n\\ \end{bmatrix}\;\; = \begin{bmatrix} 0\\ 0\\ u_i\\ v_i\\ \vdots\\ u_j\\ v_j\\ \vdots\\ u_m\\ v_m\\ \vdots\\ u_n\\ v_n\\ \vdots\\ 0\\ 0\\ \end{bmatrix}$$