有限元基础编程——杆单元（附Matlab源码）

易公子 2022年6月22日浏览：2808 评论：8 收藏：40

有限元基础编程——杆单元（附Matlab源码）

引言

”有限的单元，无限的能力“这句话来自清华大学有限元分析公开课曾攀老师的开课语。想要学好有限元这门课，不光要理解理论公式的由来及简单手酸，更要结合实际应用。本栏目将带着大家Step-By-Step基于Matlab语言实现有限元的基础操作，课程代码来自《有限元分析基础教程》——曾攀，并附赠ANSYS命令流文件进行自行验证Matlab代码正确性。

有限元“流水线套路”：

• 求解单元刚度
• 组装整体刚度
• 未知位移求解

• 本质是线性方程组求解，求解方法有很多，基于Fortran编写的可以采用JCG开源程序包，基于Matlab编写的可以采用\，默认高斯消去法，也可以使用PCG(预处理共轭梯度法)等等，总之求解线性方程组的方法很多，我们初学者可以先使用最简单的\，进行求解。
• 对于位移边界的处理，有限元有很多方法：直接消去法、置1法、拉格朗日乘子法、罚函数法等，对于初学者可以先概念最简单的直接消去法入手，等熟悉了有限元基本过程再使用更加进阶的操作。

• 节点力、应力、应变等求解

1D杆单元

这大概是最简单的有限元分析吧，简直每本有限元教材里面都会将之作为入门案例，操作虽然简单，但也是包含了有限元分析的基础步骤。

案例

有限元基础编程——杆单元（附Matlab源码）的图1

函数

% 篇幅原因，就不展示节点力、单元应力代码了，详情可在公众号(易木木响叮当)
%后台回复：杆_Matlab，即可自动获得源代码,目前百度云暂停分享功能，等分享功能
%正常时，即可上传。
function k=Bar1D2Node_Stiffness(E,A,L)
%该函数计算单元的刚度矩阵
%输入弹性模量E，横截面积A和长度L
%输出单元刚度矩阵k(2X2) 
%---------------------------------------
k=[E*A/L -E*A/L; -E*A/L E*A/L];
function z=Bar1D2Node_Assembly(KK,k,i,j)
%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k，单元的节点编号i、j
%输出整体刚度矩阵KK
%-----------------------------------
DOF(1)=i;
DOF(2)=j;
for n1=1:2
   for n2=1:2
      KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2);
   end
end
z=KK;

主程序

format long
% 典型例题[2.3(1)]P17
E1 = 2*10^5;E2 = E1;E3 = E1;
A3 = 0.06;A2 = 0.5*A3;A1 = A3/3;
L1 = 0.1;L2 = L1;L3 =L1;
k1 = Bar1D2Node_Stiffness(E1,A1,L1);
k2 = Bar1D2Node_Stiffness(E2,A2,L2);
k3 = Bar1D2Node_Stiffness(E3,A3,L3);
KK = zeros(4,4);
KK = Bar1D2Node_Assembly(KK,k1,1,2);
KK = Bar1D2Node_Assembly(KK,k2,2,3);
KK = Bar1D2Node_Assembly(KK,k3,3,4);
% 直接法求解整体刚度矩阵
k = KK([1:3],[1:3]);
p = [-100;0;50];
% 高斯消去法求解线性方程组
u = k\p

结果所得位移与手算结果一致，可自行验证。

2D杆单元

坐标变换

2D杆单元在编写的时候涉及到由局部坐标系向整体坐标系变换的过程。坐标转换矩阵 $T$ 为：

$T=\left[ \begin{matrix} \cos \alpha& \sin \alpha& 0& 0\\ 0& 0& \cos \alpha& \sin \alpha\\ \end{matrix} \right]$ 刚度矩阵、节点位移由局部坐标系 $K$ 、 $q$ 转换到整体坐标系 $\bar{K}$ 、 $\bar{q}$ ：

$\bar{K}=T^TKT$
$\bar{q}=T^Tq$ $\bar{K}=\frac{EA}{L}\left[ \begin{matrix} \cos ^2\alpha& \cos \alpha \sin \alpha& -\cos ^2\alpha& -\cos \alpha \sin \alpha\\ \cos \alpha \sin \alpha& \sin ^2\alpha& -\cos \alpha \sin \alpha& -\sin ^2\alpha\\ -\cos ^2\alpha& -\cos \alpha \sin \alpha& \cos ^2\alpha& \cos \alpha \sin \alpha\\ -\cos \alpha \sin \alpha& -\sin ^2\alpha& \cos \alpha \sin \alpha& \sin ^2\alpha\\ \end{matrix} \right]$

进行应力、节点力计算时，位移也应该由局部坐标系转换到整体坐标系中，由弹性力学中的物理方程，有1D问题的应力：

$\sigma =E\varepsilon =EBq=\frac{E}{L}\left[ \begin{matrix} -1& 1\ \end{matrix} \right] Tq$ $F=\sigma A$

案例

有限元基础编程——杆单元（附Matlab源码）的图2

函数

function k=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha)
%该函数计算单元的刚度矩阵
%输入弹性模量E，横截面积A
%输入第一个节点坐标（x1,y1），第二个节点坐标（x2,y2），角度alpha（单位是度）
%输出单元刚度矩阵k(4X4)。
L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
x=alpha*pi/180;
C=cos(x);
S=sin(x);
k=E*A/L*[C*C C*S -C*C -C*S
            C*S S*S -C*S -S*S
             -C*C -C*S C*C C*S
             -C*S -S*S C*S S*S];
function z = Bar2D2Node_Assembly(KK,k,i,j)
%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k，单元的节点编号i、j
%输出整体刚度矩阵KK
DOF(1)=2*i-1;
DOF(2)=2*i;
DOF(3)=2*j-1;
DOF(4)=2*j;
for n1=1:4
   for n2=1:4
      KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2);
   end
end
z=KK;

主程序

E=2.95e11;A=0.0001;
x1=0;y1=0;x2=0.4;y2=0;x3=0.4;y3=0.3;x4=0;y4=0.3;
alpha1=0;alpha2=90;alpha3=atan(0.75)*180/pi;
k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2,alpha1)
k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3,alpha2)  
k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3,alpha3)
k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3,alpha1)  
%建立整体刚度方程
%由于该结构共有4个节点，因此，结构总的刚度矩阵为KK(8×8)，先对K清零，
%然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装。
KK=zeros(8,8);
KK=Bar2D2Node_Assembly (KK,k1,1,2);
KK=Bar2D2Node_Assembly (KK,k2,2,3);
KK=Bar2D2Node_Assembly (KK,k3,1,3);
KK=Bar2D2Node_Assembly (KK,k4,4,3)
%边界条件的处理及刚度方程求解
k=KK([3,5,6],[3,5,6])
p=[20000;0;-25000]
u=k\p

ANSYS命令流

/ PREP7               !进入前处理  
/PLOPTS,DATE,0       !设置不显示日期和时间
!=====设置单元、材料，生成节点及单元 
ET,1,LINK1            !选择单元类型 
UIMP,1,EX, , ,2.95e11,    !给出材料的弹性模量 
R,1,1e-4,               !给出实常数(横截面积) 
N,1,0,0,0,               !生成1号节点,坐标(0,0,0)  
N,2,0.4,0,0,             !生成2号节点,坐标(0.4,0,0) 
N,3,0.4,0.3,0,            !生成3号节点,坐标(0.4,0.3,0) 
N,4,0,0.3,0,             !生成4号节点,坐标(0,0.3,0) 
E,1,2                   !生成1号单元(连接1号节点和2号节点)   
E,2,3                   !生成2号单元(连接2号节点和3号节点)
E,1,3                   !生成3号单元(连接1号节点和3号节点)
E,4,3                   !生成4号单元(连接4号节点和3号节点)
FINISH                 !前处理结束
!=====在求解模块中，施加位移约束、外力，进行求解
/SOLU                  !进入求解状态(在该状态可以施加约束及外力) 
D,1,ALL                 !将1号节点的位移全部固定 
D,2,UY，                !将2号节点的Y方向位移固定
D,4,ALL                 !将4号节点的位移全部固定 
F,2,FX,20000,             !在2号节点处施加X方向的力(20000) 
F,3,FY,-25000,             !在3号节点处施加Y方向的力(-25000) 
SOLVE                   !进行求解
FINISH                   !结束求解状态 
!=====进入一般的后处理模块 
/POST1                   !进入后处理 
PLDISP,1                  !显示变形状况
FINISH                   !结束后处理