DKT单元

前言：本文主要介绍离散基尔霍夫理论(DKT, Discrete Kirchhoff Theory)，包括其三角板单元的插值函数、刚度矩阵的推导以及薄板弯曲算例。

位移插值函数

  根据横向剪切应变板的小变形假定，板的位移场可以表述为：
$$u=-z{\phi }_x,y=-z{\phi }_y,w=w\text{\hspace{1em}(1)}$$
式中，${\phi }_x,{\phi }_y$为绕中面的转动变量，在Kirchhoff板理论中，${\phi }_x=\frac{\partial w}{\partial x},{\phi }_y=\frac{\partial w}{\partial y}$。

  如图所示，在六节点三角板单元中，节点4,5,6为各边中点。$n,\tau$分别为边界的法线和切线，$\theta$为法线$n$和$x$正向的夹角。

  对单元内转动变量进行二次插值
$${\phi }_x=\sum _{i=1}^6{N}_i{\phi }_{xi},{\phi }_y=\sum _{i=1}^6{N}_i{\phi }_{yi}\text{\hspace{1em}(2)}$$

  在各节点处，引入Kirchhoff直法线假设：
角结点处：
$$\{\begin{array}{l}{\left(\frac{\partial w}{\partial x}\right)}_i-{\phi }_{xi}=0\\ {\left(\frac{\partial w}{\partial y}\right)}_i-{\phi }_{yi}=0\end{array}(i=1,2,3)\text{\hspace{1em}(3)}$$
中点处：
$$\{\begin{array}{l}{\left(\frac{\partial w}{\partial \tau }\right)}_k-{\phi }_{\tau k}=0\\ {\phi }_{nk}=\frac{1}{2}\left({\phi }_{ni}+{\phi }_{nj}\right)\end{array}(k=4,5,6)\text{\hspace{1em}(4)}$$
  对边$ij$任意点的挠度进行$Hermite$插值得：
$$w=H_{i1}{w}_i+H_{j1}{w}_j+H_{i2}{\left(\frac{\partial w}{\partial \tau }\right)}_i+H_{j2}{\left(\frac{\partial w}{\partial \tau }\right)}_j\text{\hspace{1em}(5)}$$
式中：$H_{i1}=\left(1+\frac{2\tau }{l_{ij}}\right){\left(\frac{\tau -l_{ij}}{l_{ij}}\right)}^2$，$H_{i2}=\tau {\left(\frac{\tau -l_{ij}}{l_{ij}}\right)}^2$，$H_{j1}=\left(1-2\frac{\tau -l_{ij}}{l_{ij}}\right){\left(\frac{\tau }{l_{ij}}\right)}^2$，$H_{j2}=\left(\tau -l_{ij}\right){\left(\frac{\tau }{l_{ij}}\right)}^2$，$l_{ij}$为边$ij$的长度.。
  将式(5)两边对$\tau$求导，并将$\tau =l_{ij}/2$代入可以得到：
$${\left(\frac{\partial w}{\partial \tau }\right)}_k=-\frac{3}{2l_{ij}}{w}_i+\frac{3}{2l_{ij}}{w}_j-\frac{1}{4}{\left(\frac{\partial w}{\partial \tau }\right)}_i-\frac{1}{4}{\left(\frac{\partial w}{\partial \tau }\right)}_j\text{\hspace{1em}(6)}$$
  边界局部坐标转动变量和整体坐标转动变量存在如下关系：
$$\left[\begin{array}{c}{\phi }_n\\ {\phi }_{\tau }\end{array}\right]=\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]\left[\begin{array}{c}{\phi }_x\\ {\phi }_y\end{array}\right]\text{\hspace{1em}(7)}$$
  由式(3)、(4)、(6)以及(7)，可以得到如下关系：
$$\left[\begin{array}{cc}\cos {\theta }_{ij}& \sin {\theta }_{ij}\\ -\sin {\theta }_{ij}& \cos {\theta }_{ij}\end{array}\right]\left[\begin{array}{c}{\phi }_{xk}\\ {\phi }_{yk}\end{array}\right]=\left[\begin{array}{cccccc}0& \frac{\cos {\theta }_{ij}}{2}& \frac{\sin {\theta }_{ij}}{2}& 0& \frac{\cos {\theta }_{ij}}{2}& \frac{\sin \theta }{2}{\theta }_{ij}\\ -\frac{3}{2l_{ij}}& \frac{\sin {\theta }_{ij}}{4}& -\frac{\cos {\theta }_{ij}}{4}& \frac{3}{2l_{ij}}& \frac{\sin {\theta }_{ij}}{4}& -\frac{\cos {\theta }_{ij}}{4}\end{array}\right]{\left[\begin{array}{cccccc}{w}_i& {\phi }_{xi}& {\phi }_{yi}& w_j& {\phi }_{xj}& {\phi }_{yj}\end{array}\right]}^T(8)$$
  至此可以得到中点$k$的转动变量端点$ij$转动变量得关系：$$\left[\begin{array}{c}{\phi }_{xk}\\ {\phi }_{yk}\end{array}\right]=\left[\begin{array}{cccccc}\frac{3\sin {\theta }_{ij}}{2l_{ij}}& \frac{2{\cos }^2{\theta }_{ij}-{\sin }^2{\theta }_{ij}}{4}& \frac{3\sin {\theta }_{ij}\cos {\theta }_{ij}}{4}& -\frac{3\sin {\theta }_{ij}}{2l_{ij}}& \frac{2{\cos }^2{\theta }_{ij}-{\sin }^2{\theta }_{ij}}{4}& \frac{3\sin {\theta }_{ij}\cos {\theta }_{ij}}{4}\\ -\frac{3\cos {\theta }_{ij}}{2l_{ij}}& \frac{3\sin {\theta }_{ij}\cos {\theta }_{ij}}{4}& \frac{2{\sin }^2{\theta }_{ij}-{\cos }^2{\theta }_{ij}}{4}& \frac{3\cos {\theta }_{ij}}{2l_{ij}}& \frac{3\sin {\theta }_{ij}\cos {\theta }_{ij}}{4}& \frac{2{\sin }^2{\theta }_{ij}-{\cos }^2{\theta }_{ij}}{4}\end{array}\right]{\left[\begin{array}{cccccc}{w}_i& {\phi }_{xi}& {\phi }_{yi}& w_j& {\phi }_{xj}& {\phi }_{yj}\end{array}\right]}^T(9)$$
  根据式(9)，可以分别求得三个中点的转动变量得变换表达式，同时若令$x_{ij}=x_j-x_i,y_{ij}=y_j-y_i$，则$\sin {\theta }_{ij}=-x_{i}j/l_{ij},\cos {\theta }_{ij}=y_{ij}/l_{ij}$。所以又可以得到如下的关系：$${\left[\begin{array}{cccccc}{\phi }_{x4}& {\phi }_{y4}& {\phi }_{x5}& {\phi }_{y5}& {\phi }_{x6}& {\phi }_{y6}\end{array}\right]}^T=\mathbf{c}\mathbf{o}\mathbf{e}\mathbf{f}\bullet \mathbf{a}\text{\hspace{1em}(10)}$$
式中，$$\mathbf{a}={\left[\begin{array}{ccccccccc}{w}_1& {\phi }_{x1}& {\phi }_{y1}& w_2& {\phi }_{x2}& {\phi }_{y2}& w_3& {\phi }_{x3}& {\phi }_{y3}\end{array}\right]}^T$$
$$\mathbf{c}\mathbf{o}\mathbf{e}\mathbf{f}=\left[\begin{array}{ccccccccc}-\frac{3x_{12}}{2l_{12}^2}& \frac{2y_{12}^2-x_{12}^2}{4l_{12}^2}& -\frac{3x_{12}{y}_{12}}{4l_{12}^2}& \frac{3x_{12}}{2l_{12}^2}& \frac{2y_{12}^2-x_{12}^2}{4l_{12}^2}& -\frac{3x_{12}{y}_{12}}{4l_{12}^2}& 0& 0& 0\\ -\frac{3y_{12}}{2l_{12}^2}& -\frac{3x_{12}{y}_{12}}{4l_{12}^2}& \frac{2x_{12}^2-y_{12}^2}{4l_{12}^2}& \frac{3y_{12}}{2l_{12}^2}& -\frac{3x_{12}{y}_{12}}{4l_{12}^2}& \frac{2x_{12}^2-y_{12}^2}{4l_{12}^2}& 0& 0& 0\\ 0& 0& 0& -\frac{3x_{23}}{2l_{23}^2}& \frac{2y_{23}^2-x_{23}^2}{4l_{23}^2}& -\frac{3x_{23}{y}_{23}}{4l_{23}^2}& \frac{3x_{23}}{2l_{23}^2}& \frac{2y_{23}^2-x_{23}^2}{4l_{23}^2}& -\frac{3x_{23}{y}_{23}}{4l_{23}^2}\\ 0& 0& 0& -\frac{3y_{23}}{2l_{23}^2}& -\frac{3x_{23}{y}_{23}}{4l_{23}^2}& \frac{2x_{23}^2-y_{23}^2}{4l_{23}^2}& \frac{3x_{23}}{2l_{23}^2}& -\frac{3x_{23}{y}_{23}}{4l_{23}^2}& \frac{2x_{23}^2-y_{23}^2}{4l_{23}^2}\\ \frac{3x_{31}}{2l_{31}^2}& \frac{2y_{31}^2-x_{31}^2}{4l_{31}^2}& -\frac{3x_{31}{y}_{31}}{4l_{31}^2}& 0& 0& 0& -\frac{3x_{31}}{2l_{31}^2}& \frac{2y_{31}^2-x_{31}^2}{4l_{31}^2}& -\frac{3x_{31}{y}_{31}}{4l_{31}^2}\\ \frac{3y_{31}}{2l_{31}^2}& -\frac{3x_{31}{y}_{31}}{4l_{31}^2}& \frac{2x_{31}^2-y_{31}^2}{4l_{31}^2}& 0& 0& 0& -\frac{3y_{31}}{2l_{31}^2}& -\frac{3x_{31}{y}_{31}}{4l_{31}^2}& \frac{2x_{31}^2-y_{31}^2}{4l_{31}^2}\end{array}\right]$$
  根据式(2)和(10)得到单元内任意点的插值函数表达式为：
$${\left[\begin{array}{cc}{\phi }_x& {\phi }_y\end{array}\right]}^T=({\mathbf{N}}_1+{\mathbf{N}}_2\bullet \mathbf{c}\mathbf{o}\mathbf{e}\mathbf{f})\mathbf{a}\text{\hspace{1em}(11)}$$
式中，
$${\mathbf{N}}_1=\left[\begin{array}{ccccccccc}0& N_1& 0& 0& N_2& 0& 0& N_3& 0\\ 0& 0& N_1& 0& 0& N_2& 0& 0& N_3\end{array}\right]$$$${\mathbf{N}}_2=\left[\begin{array}{cccccc}{N}_4& 0& N_5& 0& N_6& 0\\ 0& N_4& 0& N_5& 0& N_6\end{array}\right]$$
  上述推导中，转动变量${\phi }_x,{\phi }_y$的方向和坐标轴正向并不一致，若想一致，可以参考文献[2]，这样求得的系数矩阵$\mathbf{c}\mathbf{o}\mathbf{e}\mathbf{f}$会有所不同。插值函数表达成式(11)，有利于推导和程序的实现。

单元刚度矩阵

  Kirchhoff板的几何方程为：
$${\epsilon }_{\kappa }=-z{\left[\begin{array}{ccc}\frac{{\phi }_x}{\partial x}& \frac{{\phi }_y}{\partial y}& \frac{{\phi }_x}{\partial y}+\frac{{\phi }_y}{\partial x}\end{array}\right]}^T\text{\hspace{1em}(12)}$$
  将式(10)代入几何方程，得：
$${\epsilon }_{\kappa }=-z({\mathbf{B}}_1+{\mathbf{B}}_2\bullet \mathbf{c}\mathbf{o}\mathbf{e}\mathbf{f})\mathbf{a}\text{\hspace{1em}(13)}$$
式中，
$${\mathbf{B}}_1=\left[\begin{array}{ccccccccc}0& \frac{\partial N_1}{\partial x}& 0& 0& \frac{\partial N_2}{\partial x}& 0& 0& \frac{\partial N_3}{\partial x}& 0\\ 0& 0& \frac{\partial N_1}{\partial y}& 0& 0& \frac{\partial N_2}{\partial y}& 0& 0& \frac{\partial N_3}{\partial y}\\ 0& \frac{\partial N_1}{\partial y}& \frac{\partial N_1}{\partial x}& 0& \frac{\partial N_2}{\partial y}& \frac{\partial N_2}{\partial x}& 0& \frac{\partial N_3}{\partial y}& \frac{\partial N_3}{\partial x}\end{array}\right]$$
$${\mathbf{B}}_2=\left[\begin{array}{cccccc}\frac{\partial N_4}{\partial x}& 0& \frac{\partial N_5}{\partial x}& 0& \frac{\partial N_6}{\partial x}& 0\\ 0& \frac{\partial N_4}{\partial y}& 0& \frac{\partial N_5}{\partial y}& 0& \frac{\partial N_6}{\partial y}\\ \frac{\partial N_4}{\partial y}& \frac{\partial N_4}{\partial x}& \frac{\partial N_5}{\partial y}& \frac{\partial N_5}{\partial x}& \frac{\partial N_6}{\partial y}& \frac{\partial N_6}{\partial x}\end{array}\right]$$
  单元刚度矩阵为：
$$\mathbf{K}=\int {\left({\mathbf{B}}_1+{\mathbf{B}}_2\bullet \mathbf{c}\mathbf{o}\mathbf{e}\mathbf{f}\right)}^T\mathbf{D}\left({\mathbf{B}}_1+{\mathbf{B}}_2\bullet \mathbf{c}\mathbf{o}\mathbf{e}\mathbf{f}\right)dxdy\text{\hspace{1em}(14)}$$
式中，矩阵D为：
$$\mathbf{D}=\frac{E{h}^3}{12\left(1-{\nu }^2\right)}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& \frac{1-\nu }{2}\end{array}\right]$$
  对式(14)进行换元，转换到面积坐标$({\mathbf{L}}_1,{\mathbf{L}}_2)$下求积分：
$$\mathbf{K}={\int }_0^1{\int }_0^{1-{\mathbf{L}}_2}{\left({\mathbf{B}}_1+{\mathbf{B}}_2\bullet \mathbf{c}\mathbf{o}\mathbf{e}\mathbf{f}\right)}^{\mathbf{T}}\mathbf{D}\left({\mathbf{B}}_1+{\mathbf{B}}_2\bullet \mathbf{c}\mathbf{o}\mathbf{e}\mathbf{f}\right)\left|\mathbf{J}\right|\mathbf{d}{\mathbf{L}}_1\mathbf{d}{\mathbf{L}}_2\text{\hspace{1em}(15)}$$
  整体坐标和面积坐标的关系：
$$x=\sum _{i=1}^3{L}_i{x}_i,y=\sum _{i=1}^3{L}_i{y}_i\text{\hspace{1em}(16)}$$
  所以雅克比矩阵$\mathbf{J}$为：
$$\mathbf{J}=\left[\begin{array}{cc}\frac{\partial x}{\partial L_1}& \frac{\partial y}{\partial L_1}\\ \frac{\partial x}{\partial L_2}& \frac{\partial y}{\partial L_2}\end{array}\right]=\left[\begin{array}{ccc}1& 0& -1\\ 0& 1& -1\end{array}\right]\left[\begin{array}{cc}{x}_1& y_1\\ x_2& y_2\\ x_3& y_3\end{array}\right]\text{\hspace{1em}(17)}$$
  插值函数$N_i$由面积坐标表示为：
$$\{\begin{array}{lc}{N}_i=(2L_i-1)L_i& i=1,2,3\\ N_k=4L_i{L}_j& k=4,5,6,为ij的中点\end{array}\text{\hspace{1em}(18)}$$
  根据链式求导法则，函数对于$x,y$的偏导均可以转换为对面积坐标的偏导：
$$\left[\begin{array}{c}\frac{\partial }{\partial x}\\ \frac{\partial }{\partial y}\end{array}\right]=\frac{1}{2A}\left[\begin{array}{ccc}{b}_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]{\left[\begin{array}{ccc}\frac{\partial }{\partial L_1}& \frac{\partial }{\partial L_2}& \frac{\partial }{\partial L_3}\end{array}\right]}^T\text{\hspace{1em}(19)}$$
  由式(18)和(19)就可以求得插值函数对$x,y$的偏导，至此单元刚度矩阵就可求得。关于面积坐标的内容可以查阅参考书籍[1]。

算例

  四边简支薄板受均布荷载，板尺寸为$a=b=1m,h=0.01m$，弹性模量$E=7e7$,泊松比$\nu =0.3$,均布荷载$f=1Pa$。
(1) 挠度云图

(2)中线挠度对比

  收藏点赞本文，私信博主获取算例Matlab代码。码字不易，感谢支持。
参考文献：
[1] 王勖成. 有限单元法[M]. 清华大学出版社, 2003.
[2] Batoz J L, Bathe K, Ho L W. A Study of three-node triangular bending elements[J]. International Journal for Numerical Methods in Engineering, 1980,15(12):1771-1812.