三维桁架有限元分析MATLAB代码

The finite element analysis of three-dimensional truss based on MATLAB and ABAQUS software

A 3D Truss structure has to be designed to sustain a total load 4P = 400kN, which P = 100kN. N as applied on nodes 1,2,3 and 4 as shown in Fig. 1. The foundation connecting this truss to the ground is designed support configuration I and support configuration II. Assume that the z = 0 plane represents the ground of the structure The truss members are composed of ASME A335-grade chrome steel, with the elastic modulus, yield strength, and safety factor as specified below.

E =210GPa.,SY = 205MPa.,SF = 2.5

That means the stress value in truss elements shound not exceeds 205/2.5=82MPa.

(a)	(b)

Fig. 1 Two support configurations

The support configurations I and configurations II's node numbers are shown in the following Fig. 2 and Fig. 3, respectively。

Fig. 2 The node of support configurations I

Fig. 3 The node of support configurations II

The two structures are generated using MATLAB as shown in Fig. 4.

(a)	(b)

Fig. 4 Two support configurations were implemented using MATLAB.

The MATLAB program algorithm flowchart is shown in Fig. 11.

Fig. 11 Algorithm flowchart

See matlab program for specific code.

3.1 nodal coordinates:                      

For case I:

d1=2000;d2=6000;d3=6000;h1=6000;h2=12000; %mm

%Node x-axis direction coordinates, mm

x=[-d1 d1 d1 -d1 -d2 d2 d2 -d2 -d3 d3 d3 -d3]/2;

%Node y-axis direction coordinates, mm

y=[d1 d1 -d1 -d1 d2 d2 -d2 -d2 d3 d3 -d3 -d3]/2;

%Node z-axis direction coordinates, mm

z=[h2 h2 h2 h2 h1 h1 h1 h1 0 0 0 0];

For case II:

d1=2000;d2=6000;d3=10000;h1=6000;h2=12000; %mm

%Node x-axis direction coordinates, mm

x=[-d1 d1 d1 -d1 -d2 d2 d2 -d2 -d3 d3 d3 -d3]/2;

%Node y-axis direction coordinates, mm

y=[d1 d1 -d1 -d1 d2 d2 -d2 -d2 d3 d3 -d3 -d3]/2;

%Node z-axis direction coordinates, mm

z=[h2 h2 h2 h2 h1 h1 h1 h1 0 0 0 0];

3.2 material properties:

In order to satisfy that the maximum stress in the bars does not exceed SY/SF=205/2.5=82Mpa, it is necessary to select the minimum value of the cross-section according to the stress calculated by the Matlab program. For Case 1, when the cross-section area of each rod is 1.35mm², as shown in Fig. 12 (a), the maximum stress is 81.892Mpa, which meets the requirement. For working condition 2, when the cross sectional area of each rod is 1.35mm2, as shown in Fig. 12 (b), the maximum stress is 81.892Mpa, which meets the requirement.

（a）    Case I

(b) Case II

Fig. 12 Calculation results for different structural stresses

Thus, for case I, the matlab statement is:

%Define cross-sectional area and modulus of elasticity, mm2 and MPa

A=1.35;E=2.1E5;

Thus, for case II, the matlab statement is:

A=1.35;E=2.1E5;

1.3 element connectivity matrices:

Fig. 13 Node connection number of the unit

Matlab code:

%Connection information: number, local node 1 number, local node 2 number, cross sectional area, modulus of elasticity

ele=[1 1 2 A E;2 2 3 A E;3 3 4 A E;4 1 4 A E;5 1 3 A E;6 2 4 A E;7 5 6 A E;...

  8 6 7 A E;9 7 8 A E;10 5 8 A E;11 1 5 A E;12 2 6 A E;13 3 7 A E;...

  14 4 8 A E;15 1 8 A E;16 3 8 A E;17 3 6 A E;18 1 6 A E;19 5 9 A E;...

  20 8 12 A E;21 7 11 A E;22 6 10 A E;...

  23 8 11 A E;24 7 12 A E;25 7 10 A E;26 6 11 A E;27 5 10 A E;28 6 9 A E;29 5 12 A E;30 8 9 A E];

%

3.4elemental stresses and nodal displacements:

Matlab code:

for iEle =1:EleCount

   %Calculate the displacement in local coordinates of the rod 

   n1=ele(iEle,2);n2=ele(iEle,3);

   R=CoordTransform([x(n1) x(n2)],[y(n1) y(n2)],[z(n1) z(n2)],BarLength(iEle));

   localU = R*[U(3*n1-2:3*n1,1);U(3*n2-2:3*n2,1)];

   Strain(1, iEle)=[-1/BarLength(iEle) 1/BarLength(iEle)]*localU; %Strain

   Stress(1, iEle)=ele(iEle,5)* Strain(1, iEle); %Stress

   AxialForce(1, iEle)=ele(iEle,4)* Stress(1, iEle); %AxialForce

End

%Save displacements, stresses and axial forces to a text file.

fp=fopen('Result.txt','a');

str = [char(13,10)','U',' ',num2str(U'),char(13,10)','Stress',' ',...

  num2str(Stress),char(13,10)','AxialForce',' ',num2str(AxialForce)];

fprintf(fp,str);

Therefore, for case I, the displacement, stress and axial force results are:

Therefore, for case II, the displacement, stress and axial force results are:

Fig. 17 Compared maximum stress results under different d3.