matlab的有限元编程练习

求各个节点的位移

主函数

%单元节点编号
clear all
ELEM=[1 2 5 4;
      2 3 6 5;
      4 5 8 7;
      5 6 9 8];
global C;%节点坐标矩阵
x1=-1;y1=0;
x2=0;y2=0;
x3=1;y3=0;
x4=-1;y4=1;
x5=0;y5=1;
x6=1;y6=1;
x7=-1;y7=2;
x8=0;y8=2;
x9=1;y9=2;
C=[x1,y1;x2,y2;x3,y3;x4,y4;x5,y5;x6,y6;x7,y7;x8,y8;x9,y9];
global NELEM;%单元数量
NELEM=4;
global NODE;%总节点数
NODE=9;
global NNODE;%每个单元节点数
NNODE=4;
global DOF;%节点*度
DOF=2;
%定义常量
global t;%单元厚度
t=1;%单元厚度
E=2e11;%弹性模量
v=0.2;%泊松比
global A;%单元面积
A=1;
D0=E/(1-v^2);
global D;%弹性矩阵
D=D0*[ 1 v 0;
       v 1 0;
       0 0 (1-v)/2];
%计算总体刚度
AK=AssembleK(ELEM);
%边界条件
dx1=0;dx3=0;dy3=0;dx7=0;dx9=0;dy9=0;
%载荷矩阵
syms Fx1 Fx3 Fy3 Fx7 Fx9 Fy9
F=[Fx1;0;0;0;Fx3;Fy3;0;0;0;0;0;0;Fx7;-1;0;0;Fx9;Fy9]; 
%划线法去除位移为0的总体刚度矩阵行与列
AK([1 5 6 13 17 18],:)=[];
AK(:,[1 5 6 13 17 18])=[];
%去除载荷矩阵位移为0的行
F([1 5 6 13 17 18],:)=[];
%求解位移
displacement=inv(AK)*F;
%转换为双精度浮点型
disp(['位移y1 x2 y2 x4 y4 x5 y5 x6 y6 y7 x8 y8 为'])
double(displacement)'

总体刚度矩阵函数

function AK=AssembleK(ELEM)

global NELEM;%单元数量
global NNODE;%每个单元节点数
global DOF;%节点*度
global NODE;%总节点数
global C;%节点坐标矩阵
AK=zeros(NODE*DOF,NODE*DOF);%总体刚度矩阵初始化

%总体刚度矩阵拼装
for h=1:NELEM
    hELEM=ELEM(h,:);%节点行
    EK=ELEMENTK(hELEM,C);%单元刚度矩阵
    for hh=1:NNODE
        hNUMBER=hELEM(hh);%节点总体编号
        for hhh=1:DOF
            n=DOF*(hNUMBER-1)+hhh;%总体刚度矩阵行
            c=DOF*(hh-1)+hhh;%单元刚度矩阵行
            for j=1:NNODE
                jNUMBER=hELEM(j);%节点总体编号
                for jj=1:DOF
                    m=DOF*(jNUMBER-1)+jj;%总体刚度矩阵列
                    d=DOF*(j-1)+jj;%单元刚度矩阵列
                    AK(n,m)=AK(n,m)+EK(c,d);%单元刚度矩阵加入到总体刚度矩阵
                end
            end
        end
    end
end

单元刚度矩阵函数

function EK=ELEMENTK(H,J)
global t;%单元厚度
global A;%单元面积
global D;%弹性矩阵
global NNODE;%每个单元节点数
%单元节点坐标矩阵
for k=1:NNODE
    coord(k,1)=J(H(k),1);
    coord(k,2)=J(H(k),2);
end
%定义标准单元坐标
syms a b
%插值函数
N1=(a-1)*(b-1);
co1 = subs(N1,{a,b},{-1,-1});
N1 = N1/co1;
N2=(a+1)*(b-1);
co2 = subs(N2,{a,b},{1,-1});
N2 = N2/co2;
N3=(a+1)*(b+1);
co3 = subs(N3,{a,b},{1,1});
N3 = N3/co3;
N4=(a-1)*(b+1);
co4 = subs(N4,{a,b},{-1,1});
N4 = N4/co4;
%插值函数矩阵
N = [N1,N2,N3,N4];
%插值函数的导数矩阵
for i = 1:4
	dN(1,i) = diff(N(i),a);
	dN(2,i) = diff(N(i),b);
end
%雅可比矩阵
Jacb=dN*coord;
%检验雅可比行列式是否为0
det(Jacb);
invJacb1=inv(Jacb);
%全局坐标
for p=1:NNODE
    N_xy(p,1)=invJacb1(1,:)*[dN(1,p);dN(2,p)];
    N_xy(p,2)=invJacb1(2,:)*[dN(1,p);dN(2,p)];
end
%高斯积分
aa1=1/sqrt(3);
bb1=aa1;
aa2=-aa1;
bb2=-bb1;
w1=1;%权系数
xy_GAUSS=w1*subs(N_xy,{a,b},{aa1,bb1})+w1*subs(N_xy,{a,b},{aa2,bb1})+ ...
    w1*subs(N_xy,{a,b},{aa1,bb2})+w1*subs(N_xy,{a,b},{aa2,bb2});
b1=xy_GAUSS(1,1);
b2=xy_GAUSS(2,1);
b3=xy_GAUSS(3,1);
b4=xy_GAUSS(4,1);
c1=xy_GAUSS(1,2);
c2=xy_GAUSS(2,2);
c3=xy_GAUSS(3,2);
c4=xy_GAUSS(4,2);
 %应变矩阵
B=[ b1 0 b2 0 b3 0 b4 0;
     0 c1 0 c2 0 c3 0 c4;
    c1 b1 c2 b2 c3 b3 c4 b4];
%计算单元刚度矩阵
EK=B'*D*B*t*A;

最终位移是求出来了，数量级应该没问题，但是做下来觉得单元刚度矩阵的应变矩阵求解有问题。