2D Quad4矩形单元

Q4单元介绍

平面应力问题及平面应变问题的单元列式基本一直，差别仅在弹性矩阵 $[D]$ 。

基本列式

基本方程

几何方程

对于平面应力问题，正应力 $\sigma_z=0$，且剪应力 $\tau_{zx}=\tau_{xz}=0$，$\tau_{zy}=\tau_{yz}=0$，仅考虑单元平面内的三个平面应力分量 $\sigma_x$、$\sigma_y$ 和 $\tau_{xy}=\tau_{yx}$，根据单元仅有 $\varepsilon_x$、$\varepsilon_y$ 和 $\gamma_{xy}=\gamma_{yx}$ 三个应变分量，根据弹性力学的基本公式可得平面应力问题的几何方程的矩阵形式：

$$ \left \{\varepsilon \right \} = \begin{Bmatrix}\varepsilon_x \\ \varepsilon_y \\ \varepsilon_{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}\sigma_x \\ \sigma_y \\ \tau_{xy}\end{Bmatrix}= = \begin{Bmatrix}\sigma \end{Bmatrix}=\left [ D \right ]\left \{ \varepsilon \right \}=\left [ D \right ]\begin{Bmatrix}\varepsilon_x \\ \varepsilon_y \\ \gamma_{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_x \\ \varepsilon_y \\ \gamma_{xy}\end{Bmatrix} $$

其中 $E$ 为弹性模量， $\mu$ 为泊松比。

位移场

4 节点矩形单元有 4 个节点和 4 个边，每个节点有两个平动自由度，令单元节点位移向量 $\left \{ a \right \} = \begin{Bmatrix} u_1 & v_1 & u_2 & v_2 & u_3 & v_3 & u_4 & v_4\end{Bmatrix}^T$。单元中任意一点的节点位移 $(u,v)$ 可通过线性拉格朗日插值得到：

$$ \begin{cases}u =N_1u_1+N_2u_2+N_3u_3+N_4u_4 \\ v =N_1v_1+N_2v_2+N_3v_3+N_4v_4\end{cases} $$

用矩阵表示为：

$$ \begin{Bmatrix}u \\ v\end{Bmatrix}=\begin{bmatrix}N_1 & 0 & \cdots & N_4 & 0 \\ 0 & N_1 & \cdots & 0 & N_4\end{bmatrix}\begin{Bmatrix}u_1 \\ v_1 \\ \cdots \\ u_4 \\ v_4\end{Bmatrix}=[N]\left \{ a \right \} $$

将单元节点进行编号，并用自然坐标表示，此时单元形函数 $N_i$ 为：

$$ \begin{cases}N_1(\xi,\eta) = (1-\xi)(1-\eta)/4 \\ N_2(\xi,\eta) = (1-\xi)(1-\eta)/4 \\ N_3(\xi,\eta) = (1-\xi)(1-\eta)/4 \\ N_4(\xi,\eta) = (1-\xi)(1-\eta)/4\end{cases} $$

其中， $\xi\in[-1,1]$，$\eta\in[-1,1]$。可以看出在单元节点处，$N_i=0$。

雅可比(Jacobian)矩阵

假设一般函数 $f=f(x,y)$，它是 $\xi$ 和 $\eta$ 的隐式函数（变量 $x$ 和 $y$ 均是 $\xi$ 和 $\eta$ 的函数），根据链导法则可得：

$$ \begin{cases} \frac{\partial f}{\partial \xi} = \frac{\partial f}{\partial x}\frac{\partial x}{\partial \xi} + \frac{\partial f}{\partial y}\frac{\partial y}{\partial \xi} \\ \frac{\partial f}{\partial \eta} = \frac{\partial f}{\partial x}\frac{\partial x}{\partial \eta} + \frac{\partial f}{\partial y}\frac{\partial y}{\partial \eta}\end{cases} $$

写成矩阵形式：

$$ \begin{Bmatrix}\frac{\partial f}{\partial \xi} \\ \frac{\partial f}{\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 f}{\partial x} \\ \frac{\partial f}{\partial y}\end{Bmatrix} =[J]\begin{Bmatrix}\frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y}\end{Bmatrix} $$

其中 $[J]$ 为雅可比矩阵，对于 $Q4$ 单元，可得雅可比矩阵为：

$$ [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}=\begin{bmatrix}\displaystyle \sum_{i=0}^4\frac{\partial N_i}{\partial \xi} x_i & \displaystyle \sum_{i=0}^4\frac{\partial N_i}{\partial \xi} y_i \\ \displaystyle \sum_{i=0}^4\frac{\partial N_i}{\partial \eta} x_i & \displaystyle \sum_{i=0}^4\frac{\partial N_i}{\partial \eta} y_i \end{bmatrix}=\begin{bmatrix}J_{11} & J_{12} \\ J_{21} & J_{22}\end{bmatrix} $$

则有：

$$ [J]^{-1}=\frac{1}{det([J])}\begin{bmatrix}J_{11} & -J_{12} \\ -J_{21} & J_{22}\end{bmatrix} $$

几何方程与应变矩阵

根据几何方程，将单元内任意一点的应变 $\left \{ \varepsilon \right \}$ 表示为单元节点位移 $\left \{ a \right \}$ 的函数，即：

$$ \left \{ \varepsilon \right \}=\left [ B \right ]\left \{ a \right \} $$

式中， $\left\{\varepsilon \right\}=\begin{Bmatrix}\varepsilon_x & \varepsilon_y & \tau_{xy}\end{Bmatrix}^T$，矩阵 $[B]$ 为应变矩阵。

上式还可以写成：

$$ \begin{Bmatrix}\varepsilon\end{Bmatrix}= \begin{Bmatrix}\varepsilon_x \\ \varepsilon_y \\ \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}=\frac{1}{det([J])}\begin{bmatrix}J_{22} & -J_{12} & 0 & 0 \\ 0 & 0 & -J_{21} & J_{11} \\ -J_{21} & J_{11} & J_{22} & -J_{12}\end{bmatrix}\begin{Bmatrix} \frac{\partial u}{\partial \xi} \\ \frac{\partial u}{\partial \eta} \\ \frac{\partial v}{\partial \xi} \\ \frac{\partial v}{\partial \eta} \end{Bmatrix} $$

进一步改写如下形式：

$$ \begin{Bmatrix}\varepsilon\end{Bmatrix}= \frac{1}{det([J])}\begin{bmatrix}J_{22} & -J_{12} & 0 & 0 \\ 0 & 0 & -J_{21} & J_{11} \\ -J_{21} & J_{11} & J_{22} & -J_{12}\end{bmatrix}\begin{bmatrix}\frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} & 0 \\ \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} & 0 \\ 0 & \frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} \\ 0 & \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} \end{bmatrix}\begin{Bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ u_3 \\ v_3 \\ u_4 \\v_4 \end{Bmatrix} $$

为便于计算应变矩阵 $[B]$，令 $[B]=[A][G]$，则有：

$$ [A]=\frac{1}{det([J])}\begin{bmatrix}J_{22} & -J_{12} & 0 & 0 \\ 0 & 0 & -J_{21} & J_{11} \\ -J_{21} & J_{11} & J_{22} & -J_{12}\end{bmatrix} $$

$$ [G]=\begin{bmatrix}\frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} & 0 \\ \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} & 0 \\ 0 & \frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} \\ 0 & \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} \end{bmatrix} $$

可见，对于 $Q4$ 单元，应变矩阵 $[B]$ 的元素是 $\xi$ 和 $\eta$ 的函数。

物理方程与应力矩阵

根据物理方程，单元的应力可表示为 $\{ \sigma \}=[D]\{ \varepsilon \}$。结合几何方程，可知任一点的应力可通过节点位移来表示 $\{ \sigma \} =[D][B]\{ a\}=[S]\{a\}$，其中 $[S]=[D][B]$ 为应力矩阵。

对于平面应力玩问题，弹性矩阵 $\displaystyle [D]=\frac{E}{1-u^2}\begin{bmatrix} 1 & u & 0 \\\mu & 1 & 0 \\ 0 & 0 & \frac{1-\mu}{2} \end{bmatrix}$，若为平面应变问题，弹性矩阵 $\displaystyle [D]=\frac{E(1-\mu)}{(1+\mu)(1-2\mu)}\begin{bmatrix} 1 & \frac{\mu}{1-\mu} & 0 \\\frac{\mu}{1-\mu} & 1 & 0 \\ 0 & 0 & \frac{1-2\mu}{2(1-\mu)} \end{bmatrix}$，平面问题的弹性矩阵均为 $3\times 3$ 的对称矩阵。

单元刚度矩阵

单元的应变能可表示为：

$$ U = \frac{1}{2}\int_V \left \{ \sigma \right \}^T \left \{ \varepsilon \right \} dV=\frac{1}{2}\left \{a \right \}^T \left[ \int_V [B]^T [D][B]dV \right ]\left \{ a \right \}=\frac{1}{2}\left \{ a \right \}^T\left [ t \int_{-1}^1 \int_{-1}^1 [B]^T[D][B]det([J])d\xi d\eta\right ]\left \{ a \right \} $$

由最小势能原理可得单元的刚度矩阵：

$$ [k]=t \int_{-1}^1 \int_{-1}^1 [B]^T[D][B]det([J])d\xi d\eta $$

其中 $t$ 为单元厚度。

对于 $Q4$ 单元，采用 $2\times2$ 个积分点的 Gauss 积分可求得精确结果，则刚度矩阵可表示为：

$$ [k]=t\sum_{i=1}^2\sum_{j=1}^2 w_iw_j[B(\xi_i,\eta_j)]^T[D][B(\xi_i,\eta_j)]det[J(\xi_i,\eta_j)] $$

其中 $w_i$、$w_j$ 为积分权重。

单元载荷列阵及等效节点力

单元的载荷主要包括面力（分布在单元边上的面力）、体力（分布在单位体积上的力）及温度作用等。以体力 $\{f_v\}$ 为例，给出其等效节点载荷 $\{f\}$ 的计算方法。

单元的体力表示为 $\{ f_v\} = \{ f_{vx}, f_{vy}\}^T$。单元体力 $\{ f_v\}$ 在单元变形 $\{u\}$ 上做的外力功势能可表示为：

$$ W=-\int_V\{u\}^T\{f_v\}dV $$

由最小势能原理可得，Q4 单元体力的等效节点载荷列阵为：

$$ \{f\}=\int_V[N]^T\{f_v\}dV=t\left[ \int_{-1}^1 \int_{-1}^1 [N]^Tdet([J])d\xi d\eta \right]\{f_v\} $$

与刚度矩阵类似，上式计算可通过数值积分完成，即

$$ \{ f \}=t\sum_{i=1}^2\sum_{j=1}^2w_iw_j\left [ N(\xi_i,\eta_j) \right ]^Tdet[J(\xi_i,\eta_j)]\{f_v\} $$