5 等参二维单元
5.1 概述
在前文，我们假设单元的边是直线的，即线性的，并且矩形单元的边与坐标轴平行。幸运的是，我们并不局限于只使用直线边的单元或矩形四边形单元；具有任意定向的曲线边的单元可以用来表示具有曲线边界的物体。通过使用等参单元，实现了从直线到曲线边的转换。转换过程是简单直接的；事实上，第4章中也使用了该程序，但仅限于直线边的单元，边与坐标轴平行。请记住，在使用有限元方法时必须定义两套关键关系：(1) 一套确定单元形状的关系，以及 (2) 插值函数的阶数。并不是必须在插值方程和坐标变换中使用相同的形状函数，也就是说，可以存在两套“全局”节点。在等参单元的情况下，这两套全局节点是相同的。几乎所有当前在工业和学术研究中使用的有限元程序都使用等参单元。等参单元基于具有独立变量ξ和η的局部坐标系，这些变量通常在空间坐标中从-1变化到1。相同的局部坐标关系用于定义单元内的一个变量，以及用节点的全局坐标来定义单元内某点的全局坐标。
5.2 自然坐标系
我们首先表示\(x\)和\(y\)为\(\xi\)和\(\eta\)的函数，即： \[ x=x(\xi,\eta) \] \[ \begin{equation} y=y(\xi,\eta) \end{equation} \] 其中插值函数也用\(\xi\)和\(\eta\)表示： \[\begin{equation} N_i=N_i(\xi,\eta) \end{equation} \] 由于我们在积分方程中处理的是笛卡尔导数，我们使用链式法则转换\(N_i\)的导数。 \[ \frac{\partial N_i}{\partial \xi}=\frac{\partial N_i}{\partial x}\frac{\partial x}{\partial\xi}+\frac{\partial N_i}{\partial y}\frac{\partial y}{\partial\xi}\] \[\begin{equation} \frac{\partial N_i}{\partial \eta}=\frac{\partial N_i}{\partial x}\frac{\partial x}{\partial\eta}+\frac{\partial N_i}{\partial y}\frac{\partial y}{\partial\eta} \end{equation} \] 方程3可以使用我们在前几章中使用的雅可比矩阵表达式重写，即： \[ \begin{equation} \begin{bmatrix} \frac{\partial N_i}{\partial \xi}\\\frac{\partial N_i}{\partial \eta} \end{bmatrix}=\begin{bmatrix} \frac{\partial x}{\partial\xi}&\frac{\partial y}{\partial\xi}\\ \frac{\partial x}{\partial\eta}&\frac{\partial y}{\partial\eta} \end{bmatrix} \begin{bmatrix} \frac{\partial N_i}{\partial x}\\\frac{\partial N_i}{\partial y} \end{bmatrix}=\boldsymbol{J} \begin{bmatrix} \frac{\partial N_i}{\partial x}\\\frac{\partial N_i}{\partial y} \end{bmatrix} \end{equation} \] 为了找到\(N_i\)的导数，我们将方程4转化为 \[\begin{equation} \begin{bmatrix} \frac{\partial N_i}{\partial x}\\\frac{\partial N_i}{\partial y} \end{bmatrix}=\boldsymbol{J}^{-1}\begin{bmatrix} \frac{\partial N_i}{\partial\xi}\\ \frac{\partial N_i}{\partial\eta} \end{bmatrix}=\frac{1}{|\boldsymbol{J}|}\begin{bmatrix} \frac{\partial y}{\partial\eta}&-\frac{\partial y}{\partial\xi}\\ -\frac{\partial x}{\partial\eta}&\frac{\partial x}{\partial\xi}\\ \end{bmatrix}\begin{bmatrix} \frac{\partial N_i}{\partial\xi}\\\frac{\partial N_i}{\partial\eta} \end{bmatrix} \end{equation} \] 其中\(|\boldsymbol{J}|\)表示\(\boldsymbol{J}\)的行列式。从方程5中我们得到： \[ \begin{equation} \frac{\partial N_i}{\partial x}=\frac{1}{|\boldsymbol{J}|}\left(\frac{\partial y}{\partial\eta}\frac{\partial N_i}{\partial\xi}-\frac{\partial y}{\partial\xi}\frac{\partial N_i}{\partial\eta}\right) \end{equation} \] \[ \begin{equation} \frac{\partial N_i}{\partial y}=\frac{1}{|\boldsymbol{J}|}\left(\frac{\partial x}{\partial\xi}\frac{\partial N_i}{\partial\eta}-\frac{\partial x}{\partial\eta}\frac{\partial N_i}{\partial\xi}\right) \end{equation} \] 雅可比行列式J由下式得到： \[ \begin{equation} |\boldsymbol{J}|=\frac{\partial x}{\partial\xi}\frac{\partial y}{\partial\eta}-\frac{\partial x}{\partial\eta}\frac{\partial y}{\partial\xi} \end{equation} \] 坐标系导数之间的关系为： \[ \frac{\partial\xi}{\partial x}=\frac{1}{|\boldsymbol{J}|}\frac{\partial y}{\partial\eta} \] \[ \frac{\partial\xi}{\partial y}=-\frac{1}{|\boldsymbol{J}|}\frac{\partial x}{\partial\eta} \] \[ \frac{\partial\eta}{\partial x}=-\frac{1}{|\boldsymbol{J}|}\frac{\partial y}{\partial\xi} \] \[ \begin{equation} \frac{\partial\eta}{\partial y}=\frac{1}{|\boldsymbol{J}|}\frac{\partial x}{\partial\xi} \end{equation} \] 如前所述，微分面积变换为关系： \[ \begin{equation} dA=dxdy=|\boldsymbol{J}|d\xi d\eta \end{equation} \] \(x\), \(y\)坐标由其节点值定义为： \[ x=\sum^n_{i=1}N_i[\xi,\eta]x_i \] \[ \begin{equation} y=\sum^n_{i=1}N_i[\xi,\eta]y_i \end{equation} \] 其中\(n = 4\)对于双线性四节点单元，\(n = 8\)或\(9\)对于二次或双二次四边形。假设变量\(\phi\)可以近似为 \[ \begin{equation} \phi=\alpha_1+\alpha_2x+\alpha_3y \end{equation} \] 对于节点\(i\)， \[ \begin{equation} \phi_i=\alpha_1+\alpha_2x_i+\alpha_3y_i \end{equation} \] 同样，如果我们让 \[ \begin{equation} \phi=\sum^n_{i=1}N_i\phi_i \end{equation} \] 将方程14代入方程13得到 \[ \begin{equation} \phi=\sum^n_{i=1}(\alpha_1+\alpha_2x_i+\alpha_3y_i)N_i(\xi,\eta)=\alpha_1\sum^n_{i=1}N_i(\xi,\eta)+\alpha_2\sum^n_{i=1}N_i(\xi,\eta)+\alpha_3\sum^n_{i=1}N_i(\xi,\eta) \end{equation} \] 因此 \[ \sum^n_{i=1}N_i(\xi,\eta)=1,\sum^n_{i=1}x_iN_i(\xi,\eta)=x,\sum^n_{i=1}y_iN_i(\xi,\eta)=y \]
5.3 形状函数
形状函数与第3章和第4章中用于二维(2-D)三角形和四边形的形状函数相同。 ### 5.3.1 双线性四边形 四节点双线性四边形的形状函数可以写为 \[ \begin{equation} N_i=\frac{1}{4}(1+\xi_i\xi)(1+\eta_i\eta),i=1,\cdots,4 \end{equation} \] 其中\(\xi_i\)和\(\eta_i\)代表单元边，如图1a中定义，即在节点1处，\(\xi = \eta = -1\)，因此\(\xi_01 = -1,\eta_1 = -1\); 在节点2处，\(\xi = \xi_2 = 1\), \(\eta = \eta_2 = -1\); 等等。因此展开形式方程16变成方程17： \[N_1 = (1 - \xi)(1 - \eta) / 4\] \[N_2 = (1 + \xi)(1 - \eta) / 4\] \[N_3 = (1 + \xi)(1 + \eta) / 4\] \[\begin{equation} N_4 = (1 - \xi)(1 + \eta) / 4 \end{equation}\]
在这种情况下，\(-1\le\xi\le 1\) 且\(-1\le\eta\le 1\)。如果\(\xi, \eta\)坐标系位于左下角，则形状函数变为： \[N_1 = (1 - \xi)(1 - \eta)\] \[N_2 = \xi(1 - \eta)\] \[N_3 = \xi\eta\] \[\begin{equation} N_4 = (1 - \xi)\eta \end{equation}\]
5.3.2 八节点二次四边形
对于八节点二次四边形，形状函数如下（图1b）， 角节点： \[ \begin{equation} N_i = (1 + \xi_i\xi)(1 + \eta_i\eta)(\xi_i\xi+\eta_i\eta-1) / 4 \end{equation}\] 中边节点： \[\begin{equation} \xi_i=0;N_i = (1 - \xi^2)(1 + \eta_i\eta) \end{equation}\]
5.3.3 线性三角形
将三角形单元的面积坐标系统\(（L_i, i = 1, 2, 3）\)转换为用于四边形单元的ξ, η系统（Heinrich和Pepper, 1999; Zienkiewicz, 1977; Zienkiewicz和Taylor 1989）是很简单的。三节点线性三角形的形状函数可以用\(\xi,\eta\)坐标表示为: \[ N_1=L_1=1-\xi-\eta \] \[ N_2=L_2=\xi \] \[\begin{equation} N_3=L_3=\eta \end{equation} \] 其中\(\xi,\eta\)在图2a中定义。在这种情况下\(0\xi\le 1\)且\(0\le\xi\le 1\)。
5.3.4 二次三角形
对于六节点三角形，形状函数定义为（图2b）， 角节点： \[\begin{equation}N_i = (2L_i - 1)L_i, i = 1, 2, 3\end{equation}\] 其中可以表示为\(\xi, \eta\)坐标： \[N_1=(2(1-\xi-\eta)-1)(1-\xi-\eta)\] \[N_2=(2\xi-1)\xi\] \[\begin{equation}N_3=(2\eta-1)\eta\end{equation}\] 中边节点：\[N_i = 4L_jL_k, i = 4, 5, 6; j = 2, 3, 1; k = 3, 1, 2\] 或者，以\(\xi\), \(\eta\)形式： \[N_4=4\xi\eta\] \[N_5=4\eta(1-\xi-\eta)\] \[\begin{equation}N_6=4\xi(1-\xi-\eta)\end{equation}\] 其中\(0\le\xi\le 1\)且\(0\le\eta\le 1\)。 ### 5.3.5 方向余弦 使用\(\xi\), \(\eta\)坐标，可以轻松获得与诺伊曼边界条件相关的单元边的方向余弦\(n_x\)和\(n_y\)。如图3所示的任意单元边，在点\(\boldsymbol{r}\)处的切向量\(\boldsymbol{dr}\)定义为 \[\begin{equation} \boldsymbol{dr}= dx\hat{i} + dy\hat{j} \end{equation} \] 其中\(\hat{i}\)和\(\hat{j}\)分别是\(x\)和\(y\)方向的单位向量。如果\(s\)表示沿边\(\boldsymbol{S}\)的弧长坐标，我们也可以写成 \[\begin{equation} \boldsymbol{dr}=\left(\frac{dx}{dx}\boldsymbol{\hat{i}}+\frac{dy}{ds}\boldsymbol{\hat{j}}\right)ds \end{equation} \] 其中 \[ \begin{equation} \boldsymbol{r}=x\boldsymbol{\hat{i}}+y\boldsymbol{\hat{j}} \end{equation} \] 切向量的长度为 \[ \begin{equation} |\boldsymbol{dr}|=ds\sqrt{dx^2+dy^2} \end{equation} \] 单位切向量\(\boldsymbol{\hat{t}}\)为 \[ \begin{equation} \boldsymbol{\hat{t}}=\frac{\boldsymbol{dr}}{|\boldsymbol{dr}|}=\frac{\boldsymbol{dr}}{ds}=\left(\frac{dx}{ds}\boldsymbol{\hat{i}}+\frac{dy}{ds}\boldsymbol{\hat{j}}\right) \end{equation} \] 可以通过顺时针旋转90°的旋转算子简单地获得指向单元边外侧的单位法向量。这可以通过旋转算子轻松实现 \[ \begin{equation} Rot(\theta)=\begin{bmatrix} cos\theta&-sin\theta\\ sin\theta&cos\theta \end{bmatrix} \end{equation} \] 然后，指向单元边外侧的单位法向量\(\boldsymbol{\hat{n}}\)是 \[ \begin{equation} \boldsymbol{\hat{n}}=Rot(-90\degree)\boldsymbol{\hat{t}}=\begin{bmatrix} 0&1\\-1&0 \end{bmatrix} \begin{bmatrix} \frac{dx}{ds}\\ \frac{dy}{ds} \end{bmatrix}=\frac{dy}{ds}\boldsymbol{\hat{i}}-\frac{dx}{ds}\boldsymbol{\hat{j}} \end{equation} \] 使用方程31，可以轻松获得线性和二次四边形单元以及三角形单元的方向余弦。
5.4 单元矩阵
我们将为线性四边形设置局部坐标系的矩阵关系。请记住，我们可以在单元内的任何地方设置我们的局部坐标系。一个标量变量\(\phi\)，它可以代表未知的温度、浓度或速度等，是使用方程14进行插值的。\(\phi\)的导数\(\partial\phi/\partial\xi\)和\(\partial\phi/\partial\eta\)使用下述方程表示： \[ \begin{equation} \begin{bmatrix} \frac{\partial\phi}{\partial\xi}\\ \frac{\partial\phi}{\partial\eta} \end{bmatrix}=\begin{bmatrix} \frac{\partial x}{\partial\xi}&\frac{\partial y}{\partial\xi}\\ \frac{\partial x}{\partial\eta}&\frac{\partial y}{\partial\eta} \end{bmatrix} \begin{bmatrix} \frac{\partial\phi}{\partial x}\\ \frac{\partial\phi}{\partial y} \end{bmatrix} \end{equation} \] 因此，由于\(N_i\)是用于近似\(\phi\)的形状函数，我们得到 \[ \begin{equation} \begin{bmatrix} \frac{\partial\phi}{\partial\xi}\\ \frac{\partial\phi}{\partial\eta} \end{bmatrix}=\begin{bmatrix} -(1-\eta)&(1-\eta)&\eta&-\eta\\ -(1-\xi)&-\xi&\xi&(1-\xi) \end{bmatrix} \begin{bmatrix} x_1&y_1\\ x_2&y_2\\ x_3&y_3\\ x_4&y_4\\ \end{bmatrix} \begin{bmatrix} \frac{\partial\phi}{\partial x}\\ \frac{\partial\phi}{\partial y} \end{bmatrix} \end{equation} \] 右边的前两个矩阵就是熟悉的\(2\times 2\)雅可比矩阵。
