matlab编译桁架有限元计算（附有完整代码）

最新推荐文章于 2024-02-07 14:56:09 发布

肖田变强不变秃

最新推荐文章于 2024-02-07 14:56:09 发布

阅读量9.3k

点赞数 36

文章标签：算法矩阵

本文链接：https://blog.csdn.net/hjuihui/article/details/118481823

版权

matlab编译桁架有限元计算（附有完整代码）

完整代码下载链接：
点击跳转哦下载后直接运行‘Main.m’即可

题目要求：
使用完成的代码求解的桁架示例。起重机的垂直和水平部分由铝制成（杨氏模量E=70 GPa，横截面为2 cm2）。对角桁架构件由钢制成（杨氏模量E=210 GPa，横截面为3 cm2）。如图，结构承受荷载P=6000 N。假设两个支撑节点固定（即x和y位移为0）。桁架问题
解决方案：
本问题通过使用Matlab编程实现了求解，并且运算结果与Ansys建模仿真得到的结果几乎完全一致，编程的流程如图1所示，下面我将对每一模块具体阐述。
图1 流程图

1. 对单元和节点进行编号
进行有限元分析之前，第一步就是要对单元及节点进行编号，下图2展示了我对单元和节点的具体编号，蓝色字符表示对节点的编号，红色字符代表对单元的编号，顺序为“从左到右，从上到下”
图2 节点单元编号
为了方便后续程序的读取调用，需要用将上述信息录入到Excel表格中，并且还需要包含一些关于单元和节点的其他信息，如表1中展示的是关于单元的信息，共计47个单元，其中有每个单元的起始与终端的节点编号、单元的刚度和单元的截面面积；表2中是关于节点的信息，共计25个节点，其中每个节点的横纵坐标、荷载和固定情况（表中‘1’表示固定，‘0’表示自由）

表1 单元信息
单元编号单元起始节点编号单元终端节点编号单元刚度(Pa) 单元截面面积(m2)
1 1 2 7.00E+10 2.00E-04
2 1 3 7.00E+10 2.00E-04
3 1 4 2.10E+11 3.00E-04
4 2 4 7.00E+10 2.00E-04
5 3 4 7.00E+10 2.00E-04
6 3 5 7.00E+10 2.00E-04
7 3 6 2.10E+11 3.00E-04
8 4 6 7.00E+10 2.00E-04
9 5 6 7.00E+10 2.00E-04
10 5 7 7.00E+10 2.00E-04
11 5 8 2.10E+11 3.00E-04
12 6 8 7.00E+10 2.00E-04
13 7 8 7.00E+10 2.00E-04
14 7 9 7.00E+10 2.00E-04
15 7 10 2.10E+11 3.00E-04
16 8 10 7.00E+10 2.00E-04
17 9 10 7.00E+10 2.00E-04
18 9 11 7.00E+10 2.00E-04
19 9 12 2.10E+11 3.00E-04
20 10 12 7.00E+10 2.00E-04
21 11 12 7.00E+10 2.00E-04
22 11 13 7.00E+10 2.00E-04
23 11 14 2.10E+11 3.00E-04
24 12 14 7.00E+10 2.00E-04
25 13 14 7.00E+10 2.00E-04
26 12 15 7.00E+10 2.00E-04
27 12 16 2.10E+11 3.00E-04
28 14 16 7.00E+10 2.00E-04
29 15 16 7.00E+10 2.00E-04
30 15 17 7.00E+10 2.00E-04
31 15 18 2.10E+11 3.00E-04
32 16 18 7.00E+10 2.00E-04
33 17 18 7.00E+10 2.00E-04
34 17 19 7.00E+10 2.00E-04
35 17 20 2.10E+11 3.00E-04
36 18 20 7.00E+10 2.00E-04
37 19 20 7.00E+10 2.00E-04
38 19 21 7.00E+10 2.00E-04
39 19 22 2.10E+11 3.00E-04
40 20 22 7.00E+10 2.00E-04
41 21 22 7.00E+10 2.00E-04
42 21 23 7.00E+10 2.00E-04
43 21 24 2.10E+11 3.00E-04
44 22 24 7.00E+10 2.00E-04
45 23 24 7.00E+10 2.00E-04
46 23 25 2.10E+11 3.00E-04
47 24 25 7.00E+10 2.00E-04
——————————————————————————
表2 节点信息
节点编号 x坐标(m) y坐标(m) x方向荷载(N) y方向荷载(N) x方向是否固定 y方向是否固定
1 0 0 0 0 1 1
2 1 0 0 0 1 1
3 0 2 0 0 0 0
4 1 2 0 0 0 0
5 0 4 0 0 0 0
6 1 4 0 0 0 0
7 0 6 0 0 0 0
8 1 6 0 0 0 0
9 0 8 0 0 0 0
10 1 8 0 0 0 0
11 0 10 0 0 0 0
12 1 10 0 0 0 0
13 0 11 0 0 0 0
14 1 11 0 0 0 0
15 2 10 0 0 0 0
16 2 11 0 0 0 0
17 3 10 0 0 0 0
18 3 11 0 0 0 0
19 4 10 0 0 0 0
20 4 11 0 0 0 0
21 5 10 0 0 0 0
22 5 11 0 0 0 0
23 6 10 0 0 0 0
24 6 11 0 0 0 0
25 7 11 0 -6000 0 0
————————————————————————————————
将上述的表1和表2读取到Matlab中，便于后续的计算操作，相关源码如下（对应编码文件为：input.m）：

%% 读取单元与节点数据并且统计单元和节点数量
ID_Element=xlsread('Element-node.xlsx');%读取单元编号及对应节点
ID_Node=xlsread('Node coordinates.xlsx');%读取节点坐标
Num_Element=size(ID_Element,1);%单元数量
Num_Node=size(ID_Node,1);%节点数量

2. 单元刚度矩阵的构建与全局刚度矩阵的装配
根据基础有限元知识，需要求出每一个单元对应的刚度矩阵，其计算公式为：

在这里插入图片描述
其中的每个元素都可以通过上步导入的表格数据求解，在程序中我编写了一个对应的函数来实现此功能，具体如下（对应编码文件为Element_stiffness_matrix.m）：

function [k_global,Length_element,c,s] = Element_stiffness_matrix(E,A,x1,y1,x2,y2)
%% 根据公式求解一个单元的全局刚度矩阵
Length_element=((y2-y1)^2+(x2-x1)^2)^0.5;
c=(x2-x1)/Length_element;
s=(y2-y1)/Length_element;
Transformation_matrix=[c^2 c*s -c^2 -s*c;s*c s^2 -s*c -s^2;-c^2 -s*c c^2 s*c;-s*c -s^2 s*c s^2];
k_element=A*E/Length_element;
k_global=k_element*Transformation_matrix;
end

完成此步骤之后，还需要将每个单元的刚度矩阵装配到一个总的刚度矩阵之中，其规则为根据每个单元的起始和终端编号，根据公式2，找到装配矩阵的对应位置，将其填入其中。如若存在单元矩阵重叠的情况，则让其相加。
在这里插入图片描述
对于编程而言，我采取的思路是先找到单个单元对应的位置，让其初步放到一个“装配矩阵”中，这个“装配矩阵”隶属于单个单元，所以加上了双引号。这里我也编写了一个处理此问题对应的函数，具体如下：（对应编码文件为：Assembly_stiffness_matrix.m）

function [Global_stiffness_matrix]=Assembly_stiffness_matrix(Element_stiffness_matrix,I,J,Num_node)
%% 刚度矩阵的装配 需要把单个全局刚度矩阵扩充到整体全局刚度矩阵之中 须在循环中依次进行
Global_stiffness_matrix=zeros(2*Num_node);%初始化刚度矩阵
Identifier=[2*I-1,2*I,2*J-1,2*J];
for i=1:4%通过两次循环加入到全局刚度矩阵之中
    for j=1:4
        Global_stiffness_matrix(Identifier(i),Identifier(j))=Element_stiffness_matrix(i,j);
    end
end
end

得到上述的“装配矩阵”，对每一个求和即为真正的装配矩阵，这一部分功能的实现我写在了主函数之中，对应编码如下：对应编码文件为：main.m）

clear
clc
%% 读取题目信息数据
Input
%% 单元刚度矩阵的构建与全局刚度矩阵的装配
Dimension_Matrix=2*Num_Node;%矩阵维数
K_global=zeros(Dimension_Matrix);%初始化全局矩阵
Information_element=[];%接受单元信息 长度 角度值 且按单元顺序
for i = 1:Num_Element%全局矩阵装配
    [stiffness_element,length_element,c,s]=Element_stiffness_matrix(ID_Element(i,4),ID_Element(i,5),...
        ID_Node(ID_Element(i,2),2),ID_Node(ID_Element(i,2),3),...
        ID_Node(ID_Element(i,3),2),ID_Node(ID_Element(i,3),3));
    K_global=K_global+Assembly_stiffness_matrix(stiffness_element,...
        ID_Element(i,2),ID_Element(i,3),Num_Node);
    Information_element=[Information_element;length_element,c,s];
end
%% 求解桁架整体位移以及单元固定支反力
[Displacement,Force]=Solve_problem(K_global,ID_Node(:,6),ID_Node(:,7),ID_Node(:,4),ID_Node(:,5),Num_Node);
%% 求解题目要求竖直位移
Displacement_y=[];
Displacement_x=[];
for i =1:Dimension_Matrix %竖直位移
    if mod(i,2)==0
        Displacement_y=[Displacement_y;Displacement(i)];
    else
        Displacement_x=[Displacement_x;Displacement(i)];
    end
end
[max_Displacement_y,max_Displacement_pos]=min(Displacement_y);
%% 后处理求解单元应力和应变
[Stress,Strain]=Post_treatment(Displacement,ID_Element,Information_element);
%% 求解题目最大拉压应力
[max_stress,max_pos]=max(Stress);
[min_stress,min_pos]=min(Stress);
%% 输出图像结果
Output

3. 求解桁架位移以及单元固定支反力
在本步中首先需要求解出节点位移，对于刚度矩阵需要分解一下，筛选出非固定节点对应的刚度矩阵，且外载已经给出，可以反解出非固定节点对应的位移。
在这里插入图片描述
对于编程而言，我采取的步骤是：记录非固定节点以及方向、外载组成列向量、筛选出非固定节点的对应刚度矩阵和求解出节点位移以及支反力四个步骤，同样我也对此写了一个求解问题的编程模块：（对应编程文件为：Solve_problem.m）

function [Displacement,Force]=Solve_problem(Stiffness_matrix,...
    Fixed_situation_x,Fixed_situation_y,External_load_x,External_load_y,Num_node)
%% 记录非固定节点以及方向
Record=[];%记录那一部分是固定的
for i=1:Num_node %判断是否固定x
    if Fixed_situation_x(i)==0 %只记录那些非固定的
        Record=[Record 2*i-1];
    end
    if Fixed_situation_y(i)==0
        Record=[Record 2*i];
    end
end
%% 外载组成列向量
External_load_sum=[];%外载矩阵组合成一起
for i = 1:2*Num_node
    if mod(i,2)==1
        External_load_sum(i)=External_load_x((i+1)/2);
    else
        External_load_sum(i)=External_load_y(i/2);
    end
end
%% 筛选出非固定节点的对应刚度矩阵
Free_External_load=[];%外载自由矩阵筛选
Free_stiffness=[];% 刚度矩阵筛选出自由的来
Num_record=size(Record,2);
for i =1:Num_record
    Free_External_load(i)=External_load_sum(Record(i));
    for j =1:Num_record
        Free_stiffness(i,j)=Stiffness_matrix(Record(i),Record(j));
    end
end
%% 求解出节点位移以及支反力
Free_Displacement=(Free_stiffness^-1)*Free_External_load';%求解得出自由位移
Displacement=zeros(2*Num_node,1);%整体位移 方便后续做和
for i =1:Num_record %重新组装为整体
    Displacement(Record(i))=Free_Displacement(i);
end
Force=Stiffness_matrix*Displacement;%得到各个节点所得到的支反力

end

通过上述代码，可以得到各个节点的位移以及固定节点的支反力，同时为了验证结果的准确性，我在Ansys中用相同的节点和单元编号进行的求解，对于节点位移，表3有清晰的展示：

表3 结果对比
节点编号 Matlab结果 Ansys结果
x位移(m) y位移(m) x位移(m) y位移(m)
1 0.000000 0.000000 0.000000 0.000000
2 0.000000 0.000000 0.000000 0.000000
3 0.012000 0.005143 0.012000 0.005143
4 0.012000 -0.006000 0.012000 -0.006000
5 0.046286 0.010286 0.046286 0.010286
6 0.046286 -0.012000 0.046286 -0.012000
7 0.102857 0.015429 0.102860 0.015429
8 0.102857 -0.018000 0.102860 -0.018000
9 0.181714 0.020571 0.181710 0.020571
10 0.181714 -0.024000 0.181710 -0.024000
11 0.285429 0.025714 0.285430 0.025714
12 0.282857 -0.030000 0.282860 -0.030000
13 0.345331 0.025714 0.345330 0.025714
14 0.345331 -0.032571 0.345330 -0.032571
15 0.280714 -0.095743 0.280710 -0.095743
16 0.347902 -0.095314 0.347900 -0.095314
17 0.279000 -0.165771 0.279000 -0.165770
18 0.350045 -0.165343 0.350040 -0.165340
19 0.277714 -0.239228 0.277710 -0.239230
20 0.351759 -0.238800 0.351760 -0.238800
21 0.276857 -0.315257 0.276860 -0.315260
22 0.353045 -0.314828 0.353040 -0.314830
23 0.276429 -0.393000 0.276430 -0.393000
24 0.353902 -0.392571 0.353900 -0.392570
25 0.354331 -0.471171 0.354330 0.000000
————————————————————————————————
通过上述表格对比不难看出，对同一问题不同软件的求解，得到的结果几乎一致，上述表格是经过保留6位小数的结果，得到的原数据其实有些差别，比如对于9号节点的x方向位移，Ansys的结果是0.181710，Matlab结果是0.181714，这是由于软件内部的舍入误差造成的，如果调整两者的舍入位数，可以得到完全一致的结果。
得到节点位移之后，带入全部的装配矩阵，可以求解固定节点的支反力，这个结果也可以通过理论力学的知识求得，分别对1、2节点取矩，并且列出水平和竖直平衡方程，简单计算可以得到结果。表5中列举了Matlab计算数据和手动计算得到的数据，可以看出Matlab的有些数据比如1节点x方向支反力不是0，这是由于软件内部的舍入误差造成的，而且这些数据和0十分接近，表明位移十分微小，所以这些结果的不一致并不能说明结果有误，反而更加准确的证明得到结果的有效性。

表5 固定点支反力
节点 Matlab结果手算结果
x支反力(N) y支反力(N) x支反力(N) y支反力(N)
1 -1.19545579095682e-09 -36000.0000000282 0 -36000
2 0 42000.0000000311 0 42000
——————————————————————————
4. 后处理与结果呈现
题目中还要求解单元的应力和应变，以及最大竖直位移，最大竖直位移可以很轻松得到，只要筛选出绝对值最大的位移以及对应的节点位置就可以了，对于应力和应变，我们已经得到了全局坐标系下的位移，我们要求单个单元的应力应变，即要把全局坐标系转化为局部坐标系，这需要求的每个单元对应的旋转矩阵，在第二步中，得到了每个单元的长度，正弦和余弦值，我们已经将其记录了下来，在这一步中直接应用就可以，其对应的代码我们写成了一个函数模块来求解，如下：（对应的编码文件为Post_treatment.m）

function [Stress,Strain]=Post_treatment(Displacement,ID_Element,Information_element)
%% 利用公式将全局位移转化为局部位移
Num_element=size(ID_Element,1);
Strain=[];%初始化单元局部位移
for i =1:Num_element
    u=[Information_element(i,2) Information_element(i,3) 0 0;...
        0 0 Information_element(i,2) Information_element(i,3)]...
        *[Displacement(2*ID_Element(i,2)-1);Displacement(2*ID_Element(i,2));...
        Displacement(2*ID_Element(i,3)-1);Displacement(2*ID_Element(i,3))];
    Strain=[Strain;(u(2)-u(1))/Information_element(i,1)];
end
Stress=ID_Element(:,4).*Strain;%直接运用公式求应力

end

表5 应力应变结果
单元编号 Matlab结果 Ansys结果
应力(Pa) 应变应力(Pa) 应变
1 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
2 1.8000E+08 2.5714E-03 1.8000E+08 2.5714E-03
3 8.9186E-06 4.2470E-17 4.7631E-05 2.2682E-16
4 -2.1000E+08 -3.0000E-03 -2.1000E+08 -3.0000E-03
5 -5.9501E-06 -8.5001E-17 -3.1815E-05 -4.5450E-16
6 1.8000E+08 2.5714E-03 1.8000E+08 2.5714E-03
7 9.1233E-06 4.3444E-17 4.7722E-05 2.2725E-16
8 -2.1000E+08 -3.0000E-03 -2.1000E+08 -3.0000E-03
9 -6.3144E-06 -9.0206E-17 -3.2058E-05 -4.5797E-16
10 1.8000E+08 2.5714E-03 1.8000E+08 2.5714E-03
11 1.0101E-05 4.8099E-17 4.7176E-05 2.2465E-16
12 -2.1000E+08 -3.0000E-03 -2.1000E+08 -3.0000E-03
13 -6.8001E-06 -9.7145E-17 -3.0115E-05 -4.3021E-16
14 1.8000E+08 2.5714E-03 1.8000E+08 2.5714E-03
15 1.1730E-05 5.5857E-17 4.8815E-05 2.3245E-16
16 -2.1000E+08 -3.0000E-03 -2.1000E+08 -3.0000E-03
17 -7.7716E-06 -1.1102E-16 -2.9143E-05 -4.1633E-16
18 1.8000E+08 2.5714E-03 1.8000E+08 2.5714E-03
19 1.3033E-05 6.2063E-17 4.5172E-05 2.1511E-16
20 -2.1000E+08 -3.0000E-03 -2.1000E+08 -3.0000E-03
21 -1.8000E+08 -2.5714E-03 -1.8000E+08 -2.5714E-03
22 1.4572E-06 2.0817E-17 0.0000E+00 0.0000E+00
23 1.6971E+08 8.0812E-04 1.6971E+08 8.0812E-04
24 -1.8000E+08 -2.5714E-03 -1.8000E+08 -2.5714E-03
25 3.8858E-06 5.5511E-17 0.0000E+00 0.0000E+00
26 -1.5000E+08 -2.1429E-03 -1.5000E+08 -2.1429E-03
27 -2.8284E+07 -1.3469E-04 -2.8284E+07 -1.3469E-04
28 1.8000E+08 2.5714E-03 1.8000E+08 2.5714E-03
29 3.0000E+07 4.2857E-04 3.0000E+07 4.2857E-04
30 -1.2000E+08 -1.7143E-03 -1.2000E+08 -1.7143E-03
31 -2.8284E+07 -1.3469E-04 -2.8284E+07 -1.3469E-04
32 1.5000E+08 2.1429E-03 1.5000E+08 2.1429E-03
33 3.0000E+07 4.2857E-04 3.0000E+07 4.2857E-04
34 -9.0000E+07 -1.2857E-03 -9.0000E+07 -1.2857E-03
35 -2.8284E+07 -1.3469E-04 -2.8284E+07 -1.3469E-04
36 1.2000E+08 1.7143E-03 1.2000E+08 1.7143E-03
37 3.0000E+07 4.2857E-04 3.0000E+07 4.2857E-04
38 -6.0000E+07 -8.5714E-04 -6.0000E+07 -8.5714E-04
39 -2.8284E+07 -1.3469E-04 -2.8284E+07 -1.3469E-04
40 9.0000E+07 1.2857E-03 9.0000E+07 1.2857E-03
41 3.0000E+07 4.2857E-04 3.0000E+07 4.2857E-04
42 -3.0000E+07 -4.2857E-04 -3.0000E+07 -4.2857E-04
43 -2.8284E+07 -1.3469E-04 -2.8284E+07 -1.3469E-04
44 6.0000E+07 8.5714E-04 6.0000E+07 8.5714E-04
45 3.0000E+07 4.2857E-04 3.0000E+07 4.2857E-04
46 -2.8284E+07 -1.3469E-04 -2.8284E+07 -1.3469E-04
47 3.0000E+07 4.2857E-04 3.0000E+07 4.2857E-04
——————————————————————————————
为了更加直观的展示结果，我们以真实图像的方式来展示，下图分别是位移偏离结果，最大拉压应力以及最大竖直位移的图像，其中位移偏离图也有Ansys的仿真结果，至此，我们已经完全求解出各种数据以及图像结果。
matlab计算位移图
ansys计算位移图
matlab计算最大应力
matlab计算最大垂直位移
从计算速度上来看，Matlab代码运行时间较长，大概十秒左右，而Ansys几乎瞬间求解完成，我自己编译的Matlab代码效率不是很高，冗杂性造成了所占计算机内存过大，具有很大的优化空间。

肖田变强不变秃

关注

36
点赞
踩
134

收藏

觉得还不错? 一键收藏
15
评论
matlab编译桁架有限元计算（附有完整代码）

matlab编译桁架有限元计算（附有完整代码）**题目要求：**使用完成的代码求解的桁架示例。起重机的垂直和水平部分由铝制成（杨氏模量E=70 GPa，横截面为2 cm2）。对角桁架构件由钢制成（杨氏模量E=210 GPa，横截面为3 cm2）。如图，结构承受荷载P=6000 N。假设两个支撑节点固定（即x和y位移为0）。解决方案本问题通过使用Matlab编程实现了求解，并且运算结果与Ansys建模仿真得到的结果几乎完全一致，编程的流程如图1所示，下面我将对每一模块具体阐述。1. 对单元和节点进行编
复制链接

扫一扫