有限元基础编程——Q4单元

易公子 2022年7月23日浏览：2267 评论：4 收藏：9

有限元基础编程——Q4单元

本次分享的是2D四节点矩形单元有限元编程，首先介绍等参单元的概念，然后围绕着编程相关的理论进行讲解，结合简单案例编制相应的Matlab程序进行求解，与Abaqus求解结果(位移)进行对比，验证程序的正确性。程序中处理边界条件的方法是“置1法”，将数据进行预处理，更加适应Abaqus或者ANSYS等CAE仿真软件前处理得到输入数据文件。代码部分风格借鉴了崔济东老师的《有限单元法——编程与应用》，运用两点高斯积分对单元刚度进行求解，避免了曾攀老师书中的Matlab符号解（耗时较大）。程序、Abaqus计算文件等全部内容可在公众号：易木木响叮当，回复：“Q4单元Matlab” 即可自动获取。

平面四节点单元等参描述

等参单元的概念相信大家已经多多少少从有限元教材中了解了一点点，木木再带着大家重温一遍，讲述其容易迷惑的地方。

有限元基础编程——Q4单元的图1

我们可以从上图看出，任意四边形被边长为 1 的正方形所代替，坐标系也相应地发生了改变，改变后的坐标系习惯上称为“母单元坐标系”，母单元的形函数我们是已知的（前辈们推导的），那么我们在参与计算时还需要知道原来坐标系“子单元坐标系”的形函数，这时候出现了一个重要的概念——雅可比(Jacobian)矩阵，本身是一个数学概念，用于有限元中表示两个坐标系的转换，即：

$\left[ \begin{array}{c} \frac{\partial N}{\partial \xi}\\ \frac{\partial N}{\partial \eta}\\ \end{array} \right] =J\left[ \begin{array}{c} \frac{\partial N}{\partial x}\\ \frac{\partial N}{\partial y}\\ \end{array} \right]$

对于平面四节点单元，Jacobian矩阵可写为：

$\left[ J \right] =\left[ \begin{matrix} \frac{\partial x}{\partial \xi}& \frac{\partial y}{\partial \xi}\\ \frac{\partial x}{\partial \eta}& \frac{\partial y}{\partial \eta}\\ \end{matrix} \right] =\left[ \begin{matrix} \sum_{i=1}^4{\frac{\partial N_i}{\partial \xi}x_i}& \sum_{i=1}^4{\frac{\partial N_i}{\partial \xi}y_i}\\ \sum_{i=1}^4{\frac{\partial N_i}{\partial \eta}x_i}& \sum_{i=1}^4{\frac{\partial N_i}{\partial \eta}y_i}\\ \end{matrix} \right] =\left[ \begin{matrix} J_{11}& J_{12}\\ J_{21}& J_{22}\\ \end{matrix} \right]$ 从以上两个式子对Jacobian矩阵的描述，我们可以看出：Jocabian矩阵本身是表示两个坐标系Mapping的关系；可以通过Jacobian进行形函数偏导数的转换。

单元刚度求解

本次分享只讲述计算模型未知位移（基于位移的有限元操作核心）求解，其余未知量可通过线弹性转换关系进行计算。

应变矩阵 B

单元任一点应变可通过应变矩阵 B与单元结点位移相联系：

$\varepsilon =Bq^e$ 为了便于编程上式可写为：

$\varepsilon =\left[ \begin{array}{c} \varepsilon _x\\ \varepsilon _y\\ \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\\ \end{array} \right] =\frac{1}{\det J}\left[ \begin{matrix} J_{22}& -J_{12}& 0& 0\\ 0& 0& -J_{21}& J_{11}\\ -J_{21}& J_{11}& J_{22}& -J_{12}\\ \end{matrix} \right] \left[ \begin{matrix} \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{matrix} \right] \left[ \begin{array}{c} u_1\\ v_1\\ u_2\\ v_2\\ u_3\\ v_3\\ u_4\\ v_4\\ \end{array} \right]$

两点高斯——勒让德积分

有限元基础编程——Q4单元的图2

从上图我们可以看出单元积分点的分布特别符合高斯——勒让德积分域（-1，1），于是采用高斯——勒让德积分，相应的积分点，权系数数值分析教材中肯定会涉及，针对两点高斯积分，积分点为 $\pm \frac{1}{\sqrt{3}}$ ，权系数为1。

$k=t\sum_{i=1}^2{\sum_{j=1}^2{w_iw_jB^TDB}\det \left[ J \right]}$

模型实例

有限元基础编程——Q4单元的图3

模型还是采用上次分享的CST单元编程使用的矩形薄板，此节点顺序为曾攀老师书中的节点顺序，比较适用“降阶法”，此次程序的数据文件由Abaqus产生，节点顺序也相应由Abaqus产生。

有限元基础编程——Q4单元的图4

数据准备

输入信息文件可由Abaqus产生，也可由任意具有前处理功能的软件产生，因木木个人喜好Abaqus, 故本文数据由 Abaqus 产生，可由Inp文件Copy节点、单元信息。

% 材料参数
E = 1e7;
NU = 1/3;
t = 0.1;
ID = 1;
% 输入信息文件
node = [0 0
          1 0
          2 0
          0 1
          1 1
          2 1];
element = [1, 2, 5, 4
              2, 3, 6, 5];
force = [3 2 -5000
            6 2 -5000];
constrain = [1 1 0
                 1 2 0
                 4 1 0
                 4 2 0];
%确定节点、单元个数
[nnode,ntmp]=size(node);
[nelem,etmp]=size(element);
[nforce,ftmp]=size(force);
[nconstrain,ctmp]=size(constrain);

%预先设定总体刚度矩阵、节点力向量、节点约向量
%针对二维问题
KKG=zeros(2*nnode);
FFG=zeros(2*nnode,1);
UUG=zeros(2*nnode,1);
StressElem=zeros(nelem,3);
StressNode=zeros(nnode,3);
k=zeros(8,8);     % 单元刚度矩阵

置 1 法

本次处理边界条件时，摒弃了以往的“划行划列”式的计算，每次计算位移都需要调整刚度矩阵的维度，很不方便。

原理

设未知位移自由度处 r 对应的整体刚度矩阵 $K_{rr}$ 为 1 ，同时也置对应的载荷元素 $p_r=0$ ，所有未知位移进行循环，这样得出的刚度方程右边的载荷项就是已知的，直接对刚度方程进行求解即可。具体可看《有限元基础教程》P178.

相应程序

% 置1法处理边界条件
for i=1:nforce
    m=force(i,1);
    n=force(i,2);
    P(2*(m-1)+n)=force(i,3);
end

for i = 1:nconstrain
    m=constrain(i,1);
    n=constrain(i,2);
    U(2*(m-1)+n)=constrain(i,3);
    KKG(2*(m-1)+n,:)=0;
    KKG(:,2*(m-1)+n)=0;
    KKG(2*(m-1)+n,2*(m-1)+n)=1;
    P(2*(m-1)+n)=0;
end
U=KKG\P

Abaqus计算验证

有限元基础编程——Q4单元的图5

节点序号	U1		U2
节点序号	Matlab	Abaqus	Matlab	Abaqus
1	0	-20.E-33	0	-5.E-33
2	-0.0600	-0.0600	-0.0867	-0.0867
3	-0.0800	‍-0.0800	-0.2533	-0.0867
4	0	-20.E-33	0	-5.E-33
5	0.0600	0.0600	-0.0867	-0.0867
6	0.0800	0.0800	-0.2533	-0.2533