有限元的学习笔记

Sink Arsenic

已于 2022-12-24 21:15:21 修改

阅读量8.3k

点赞数 10

分类专栏：学习笔记文章标签：学习 matlab 有限元

于 2022-11-25 22:28:03 首次发布

本文链接：https://blog.csdn.net/qq_44846756/article/details/128045652

版权

学习笔记专栏收录该内容

22 篇文章 5 订阅

订阅专栏

杂七杂八

知识点
 有限元刚度矩阵推导
 例题
 ppt
matlab 特征值特征向量计算
 西工大题
 轴对称问题有限元法
 高斯-赛德尔（Gauss-Seidel）迭代法求解线性方程组

A=(a1,a2)
B=(b1,b2)

A×B
=(a1,a2)×(b1,b2)
=a1b2-a2b1

A=(a1,a2,a3)
B=(b1,b2,b3)

A×B
=(a1,a2,a3)×(b1,b2,b3)
=(a2b3-a3b2,a3b1-a1b3,a1b2-a2b1)

什么是应力？

一个物体在外力作用下要产生变形，如果要将这个物体分成几个部分，那么被分开部分的表面之间就有相互作用力来反抗物体的变形，这种物体内部的作用力为内力，也就是我们通常所说的应力。

什么是面力和体力？

一个物体所受到的外力可分为两种：一种是分布在物体表面的作用力，称为面力，例如一个物体对另一个物体的压力，水压力等。另一种是分布在物体体积内的力，称为体力，如重力，磁场力，惯性力等。

什么是应力矢量？

从一个处于平衡状态的的物体中分离出一小部分B，这小部分的面积为S，在S上任取一点Q，在Q点的邻域

△S上作用有合力△Fv，其中v是小部分B的过Q点切平面的外法线，那么当
在这里插入图片描述

即应力矢量应力分量应力张量

那么应力矢量应力分量应力张量
就称为应力矢量。

要注意的是：1.应力矢量表示位于B部分S面上Q点，且其切平面的法向为v的应力矢量，即使Q点不动，但是v的方向改变时所切的物体B也改变，则应力矢量也要改变。

应力矢量是一个固定矢量

3.应力矢量可在过Q点的切平面上进行分解成剪应力分量和正应力分量。

什么是应力分量？

从前面的讲解可以看出，应力分量就是应力矢量的分量。如果选取一个笛卡尔坐标面为过Q点的一个切平面，其上的应力矢量可以沿坐标方向进行分解成一个正应力分量和两个剪应力分量。这些应力分量的方向定义如下：当所选的坐标面的法向方向与坐标轴的正向相同时，规定各个应力分量的正向与坐标轴的指向相同，反之，当所选的坐标面的法向方向指向坐标轴的负向时，规定各个应力分量的正向也指向坐标轴的负向。

什么是应力张量？

在笛卡尔坐标系中过Q点可以做相互正交的三个坐标面，那么这三个坐标面上总共就有九个应力分量，它们是一个整体，它们共同构成了一个新的物理量，即就是应力张量应力矢量应力分量应力张量
，其中i表示应力分量所在的坐标面的外法线的指向，j表示应力分量的方向。i和j通常表示x，y，z，但经常用1，2，3 来表示，这样，这九个应力分量就可以表示成一个矩阵为：

应力矢量应力分量应力张量

平面问题有限元分析四节点矩形单元.
边界约束的处理

有限元题目 and 代码

10
在这里插入图片描述

syms a b x y;
A(a,b)=[1/8*(-6*x/a^3)*(1+y/b),0,a/8*(1+y/b)*(-2/a^2-6*x/a^3);1/8*(-6*y/b^3)*(1+x/a),-b/8*(1+x/a)*(-2/b^2-6*y/b^3),0;1/4*(2/(a*b)+1/b*(1/a-3*x^2/a^3)+1/a*(1/b-3*y^2/b^3)),-b/(4*a)*(1/b-2*y/b^2-3*y^2/b^3),a/(4*b)*(1/a-2*x/a^2-3*x^2/a^3)];
B=[A(-1/2,-1/2) A(1/2,-1/2) A(1/2,1/2) A(-1/2,1/2)];
E=200*10^9;
u=0.3;
D=E/(1-u^2)*[1,u,0;u,1,0;0,0,(1-u)/2];
t=0.01;
k=t^3/12*B.'*D*B;
ky=int(k,y,-1/2,1/2);
K=int(ky,x,-1/2,1/2);
Ke=zeros(27,27);
%贡献自由度
%分离
for i=1:3
    for j=1:3
        G11(i,j)=K(i,j);
        G12(i,j)=K(i,j+3);
        G13(i,j)=K(i,j+6);
        G14(i,j)=K(i,j+9);
        G21(i,j)=K(i+3,j);
        G31(i,j)=K(i+6,j);
        G41(i,j)=K(i+9,j);
        G22(i,j)=K(i+3,j+3);
        G23(i,j)=K(i+3,j+6);
        G24(i,j)=K(i+3,j+9);
        G32(i,j)=K(i+6,j+3);
        G42(i,j)=K(i+9,j+3);
        G33(i,j)=K(i+6,j+6);
        G34(i,j)=K(i+6,j+9);
        G43(i,j)=K(i+9,j+6);
        G44(i,j)=K(i+9,j+9);
    end
end
for i=1:3
    for j=1:3
        %1
        Ke(i,j)=G11(i,j)%11
        Ke(i,j+3)=G12(i,j)%12
        Ke(i,j+9)=G14(i,j)%14
        Ke(i,j+12)=G13(i,j)%15
        Ke(i+3,j)=G21(i,j)%21
        Ke(i+9,j)=G41(i,j)%41
        Ke(i+12,j)=G31(i,j)%51
        %2
        Ke(i+3,j+3)=G22(i,j)+G11(i,j)%22
        Ke(i+3,j+6)=G12(i,j)%23
        Ke(i+6,j+3)=G21(i,j)%32
        Ke(i+3,j+9)=G24(i,j)%24
        Ke(i+9,j+3)=G42(i,j)%42
        Ke(i+3,j+12)=G23(i,j)+G14(i,j)%25
        Ke(i+12,j+3)=G32(i,j)+G41(i,j)%52
        Ke(i+3,j+15)=G13(i,j)%26
        Ke(i+15,j+3)=G31(i,j)%26
        %3
        Ke(i+6,j+6)=G22(i,j)%33
        Ke(i+6,j+12)=G24(i,j)%35
        Ke(i+12,j+6)=G42(i,j)%53
        Ke(i+6,j+15)=G23(i,j)%36
        Ke(i+15,j+6)=G32(i,j)%63
        %4
        Ke(i+9,j+9)=G44(i,j)+G11(i,j)%44
        Ke(i+9,j+12)=G43(i,j)+G12(i,j)%45
        Ke(i+12,j+9)=G34(i,j)+G21(i,j)%54
        Ke(i+9,j+18)=G14(i,j)%47
        Ke(i+18,j+9)=G41(i,j)%74
        Ke(i+9,j+21)=G13(i,j)%48
        Ke(i+21,j+9)=G31(i,j)%84
        %5
        Ke(i+12,j+12)=G11(i,j)+G22(i,j)+G33(i,j)+G44(i,j)%55
        Ke(i+12,j+15)=G43(i,j)+G12(i,j)%56
        Ke(i+15,j+12)=G34(i,j)+G21(i,j)%56
        Ke(i+12,j+18)=G24(i,j)%57
        Ke(i+18,j+12)=G42(i,j)%75
        Ke(i+12,j+21)=G14(i,j)+G23(i,j)%58
        Ke(i+21,j+12)=G41(i,j)+G32(i,j)%85
        Ke(i+12,j+24)=G13(i,j)%59
        Ke(i+24,j+12)=G31(i,j)%95
        %6
        Ke(i+15,j+15)=G22(i,j)+G33(i,j)%66
        Ke(i+15,j+21)=G24(i,j)%68
        Ke(i+21,j+15)=G42(i,j)%86
        Ke(i+15,j+24)=G23(i,j)%69
        Ke(i+24,j+15)=G32(i,j)%69
        %7
        Ke(i+18,j+18)=G44(i,j)%77
        Ke(i+18,j+21)=G43(i,j)%78
        Ke(i+21,j+18)=G34(i,j)%87
        %8
        Ke(i+21,j+21)=G44(i,j)+G33(i,j)%88
        Ke(i+21,j+24)=G43(i,j)%89
        Ke(i+24,j+21)=G43(i,j)%98
        %9
        Ke(i+24,j+24)=G33(i,j)%99
    end
end

9
在这里插入图片描述

E= 200*10^9;
A= 9*10^-3;
L1 = 2;
L2 = 1;
iz = 3*10^-5;
a1 = -pi/2;
a2 = 0;
P = 10000;

% E= 210*10^9;
% A= 9*10^-3;
% L = 4;
% iz = 7.5*10^-6;
% a = pi/2;

b1 = E*A/L1;
b2 = E*A/L2;
c1 = 12*E*iz/L1^3;
d1 = 6*E*iz/L1^2;
e1 = 4*E*iz/L1;
c2 = 12*E*iz/L2^3;
d2 = 6*E*iz/L2^2;
e2 = 4*E*iz/L2;


k1 = [b1 0 0 -b1 0 0; 0 c1 d1 0 -c1 d1; 0 d1 e1 0 -d1 e1/2;
        -b1 0 0 b1 0 0; 0 -c1 -d1 0 c1 -d1; 0 d1 e1/2 0 -d1 e1]

k2 = [b2 0 0 -b2 0 0; 0 c2 d2 0 -c2 d2; 0 d2 e2 0 -d2 e2/2;
        -b2 0 0 b2 0 0; 0 -c2 -d2 0 c2 -d2; 0 d2 e2/2 0 -d2 e2]

T1 = [cos(a1) -sin(a1) 0 0 0 0; sin(a1) cos(a1) 0 0 0 0;
       0 0 1 0 0 0; 0 0 0 cos(a1) -sin(a1) 0; 
       0 0 0 sin(a1) cos(a1) 0; 0 0 0 0 0 1]
   
T2 = [cos(a2) -sin(a2) 0 0 0 0; sin(a2) cos(a2) 0 0 0 0;
       0 0 1 0 0 0; 0 0 0 cos(a2) -sin(a2) 0; 
       0 0 0 sin(a2) cos(a2) 0; 0 0 0 0 0 1]

K1 = T1*k1*T1'
K2 = T2*k2*T2'
K3 = K2
K = zeros(12,12);
for i = 1:6
    for j = 1:6
        K(i,j) = K1(i,j);
    end
end
for i = 1:6
    for j = 1:6
        K(i+3,j+3) = K(i+3,j+3)+K2(i,j);
    end
end
for i = 1:6
    for j = 1:6
        K(i+6,j+6) = K(i+6,j+6)+K3(i,j);
    end
end

Ks = zeros(9,9);
for i = 1:9
    for j = 1:9
        Ks(i,j) = K(i,j);
    end
end


Kb = [-P;0;0;0;0;0;0;-2*P;0];
x = Ks\Kb;
x(10:12)=0;
fxym = zeros(12,1);
for i = 1:12
     fxym(i) = K(i,1:12)*x(1:12)
end

syms x 
l = 2;
N1 = 1-x/l;
N2 = x/l;
N3 = 1-3/l^2*x^2+2/l^3*x^3;
N4 = x-2/l*x^2+1/l^2*x^3;
N5 = 3/l^2*x^2-2/l^3*x^3;
N6 = -1/l*x^2+1/l^2*x^3;

w = [-0.0211;-0.0094;0.01167;
         -1.11e-05;-0.0094;0.0083;-5.5556e-06;-0.0028;0.0050;0;0;0];

[N1 0 0 N2 0 0; 0 N3 N4 0 N5 N6]*w(1:6)
a = 5.999999*10^6;
b = 3.33*10^-7;
c = 2*10^11;
d1 = 0.0105;
d2 = 5.54*10^-6;
d3 = 5.56*10^-6;

x=linspace(0,1);
% y1 = (0.03*x-0.0158)*a*b+c*d1;
% plot(x,y1)

% y2 = (6*10^-4*x-0.0018)*a*b+c*d2;
% plot(x,y2)
y3 = (-36*10^-4*x-0.0016)*a*b+c*d3;
plot(x,y3)