有限元基础编程——杆单元(附Matlab源码)
2022-07-28 本文已影响0人
易木木响叮当
本篇内容转载本人公众号:易木木响叮当,涉及代码可以在后台回复:杆_Matlab,即可自动获取。
引言
”有限的单元,无限的能力“这句话来自清华大学有限元分析公开课曾攀老师的开课语。想要学好有限元这门课,不光要理解理论公式的由来及简单手酸,更要结合实际应用。本栏目将带着大家Step-By-Step基于Matlab语言实现有限元的基础操作,课程代码来自《有限元分析基础教程》——曾攀,并附赠ANSYS命令流文件进行验证Matlab代码正确性。
有限元“流水线套路”:
-
求解单元刚度
-
组装整体刚度
-
未知位移求解
-
本质是线性方程组求解,求解方法有很多,基于Fortran编写的可以采用JCG开源程序包,基于Matlab编写的可以采用
\
,默认高斯消去法,也可以使用PCG(预处理共轭梯度法)等等,总之求解线性方程组的方法很多,我们初学者可以先使用最简单的\
,进行求解。 -
对于位移边界的处理,有限元有很多方法:直接消去法、置1法、拉格朗日乘子法、罚函数法等,对于初学者可以先概念最简单的直接消去法入手,等熟悉了有限元基本过程再使用更加进阶的操作。
-
-
节点力、应力、应变等求解
1D杆单元
这大概是最简单的有限元分析吧,简直每本有限元教材里面都会将之作为入门案例,操作虽然简单,但也是包含了有限元分析的基础步骤。
案例
例题.jpg函数
function k=Bar1D2Node_Stiffness(E,A,L)
k=[E*A/L -E*A/L; -E*A/L E*A/L];
function z=Bar1D2Node_Assembly(KK,k,i,j)
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杆单元在编写的时候涉及到由局部坐标系向整体坐标系变换的过程。坐标转换矩阵为:
刚度矩阵、节点位移由局部坐标系、转换到整体坐标系、:
进行应力、节点力计算时,位移也应该由局部坐标系转换到整体坐标系中,由弹性力学中的物理方程,有1D问题的应力:
案例
杆单元案例函数
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 !结束后处理