MATLAB CODE FOR 3D TRUSS
A matlab code is developed to solve the 3D truss problem as given in D.A Hutton Fundamental of finite element analysis. The problem is analyzed in ansys to verify the results. Problem The three member truss is connected by a ball and socket joint and fixed at nodes 1,3 and 4 as shown in the figure below. A 5000-lb is force is applied in negative Y direction at node 4.The area, length and modulus of elasticity of each element is 0.5 in 2,50 in and 30x106 psi. Compute the displacements at node 4.
Fig (a) shows the coordinates of nodes in global coordinate system and the load acting on the structure. Fig(b) Shows the nodal and element numbering Matlab code clear all; clc;
% NCA(Nodal collection array) is 2d array consisting of element number ,node % numbers,area,Young's modulus, modulus, coordinates of nodes corresponding to the % each element in row wise.
1
NCA=[1 1 4 0.5 30*10^6 50 0 0 30 40 0 0;2 2 4 0.5 30*10^6 50 0 0 -30 40 0 0 ;3 3 4 0.5 30*10^6 50 0 -30 0 40 0 0]; for EN=1:3 A=NCA(EN,4); EM=NCA(EN,5); L=NCA(EN,6); K=A*EM/L; % LSTIFF is stiffness matrix(2x2) element LSTIFF=[K -K;-K K]; % TM is transformation matrix of elements TM(:,:,EN) =[(NCA(EN,10)-NCA(EN,7))/NCA(EN,6) (NCA(EN,11)- NCA(EN,8))/NCA(EN,6) (NCA(EN,12)-NCA(EN,9))/NCA(EN,6) 0 0 0;0 0 0 (NCA(EN,10)-NCA(EN,7))/NCA(EN,6) (NCA(EN,11)-NCA(EN,8))/NCA(EN,6) (NCA(EN,12)-NCA(EN,9))/NCA(EN,6)]; % ESTIFF is stiffness matrix stored as 3D matrix ESTIFF(:,:, EN)=(TM(:,:,EN))'*LSTIFF*(TM(:,:,EN)); % Row numbering and column numbering relating element stiffness matrix and global stiffness matrix through nodal numbering nodal numbering for NN=1:2 RN(EN,3*NN-2)=3*NCA(EN,NN+1)-2; RN(EN,3*NN-1)=3*NCA(EN,NN+1)-1; RN(EN,3*NN)=3*NCA(EN,NN+1); end for NN=1:2 CN(EN,3*NN-2)=3*NCA(EN,NN+1)-2; CN(EN,3*NN-1)=3*NCA(EN,NN+1)-1; CN(EN,3*NN)=3*NCA(EN,NN+1); end end % GSTIFF Global stiffness matrix % GSTIFF initialized to zero GSTIFF=zeros(12); for EN=1:3 for I=1:6 for J=1:6 GSTIFF(RN(EN,I),CN(EN,J))=GSTIFF(RN(EN,I),CN(EN,J))+ESTIFF(I,J,EN); end end end % boundary conditions % The constrained degree of freedom is assigned as 1 and free degree of freedom as 0 2
stat=[1 1 1 1 1 1 1 1 1 0 0 0]; % KK – Stiffness matrix after imposing boundary conditions p=0; for i=1:12 if(~stat(i)) p=p+1; q=0; for j=1:12 if(~stat(j)) q=q+1; KK(p,q)=GSTIFF(i,j); end end end end % Force acting on the truss corresponding to displacements to be determine FF=[0 -5000 0]'; % solving the equation for displacements disp('displacements'); qq=KK\FF; disp(qq); % stress and forces calculation E=30*10^6; for EN=1:3 % Stresses in elements stress(:,EN)=E*[-1/NCA(EN,6) 1/NCA(EN,6)]*TM(:,:,EN)*[0 0 0 qq(1) qq(2) qq(3)]'; % Forces at nodes Force(:,EN)=ESTIFF(:,:,EN)*[0 0 0 qq(1) qq(2) qq(3)]'; % Resultant force at nodes force(:,EN)= TM(:,:,EN)*ESTIFF(:,:,EN)*[0 0 0 qq(1) qq(2) qq(12)]'; end EN=1:3; disp('Stresses'); fprintf('%f\n\n',stress(:,EN)); disp('Resultant forces at nodes'); fprintf('%f\n\n',force(:,EN)); disp('Forces at nodes'); fprintf('%f\n\n',Force(:,EN));
3
ANSYS procedure
The key steps involved in analysis the model using the ANSYS is detailed below Preprocessing 1.Specify title 2.Specify preferences 3.Define elements 4.Define real contants 5.Define material properties 6.Define nodes 7. Define elements Solution 8. Apply displacements 9. Apply loads 10.Solve Post processing 11.Read results 12. Print deformed results 13. Exit ansys Preprocessing 1.Specify title Utility menu: File Change title a. Enter the title ‘Static analysis of 3D truss’ b. Click on ok 2. Specify preferences Main menu: Preferences a. Turn on Structural b. Click on ok 3. Define element type Main menu: Preprocessor Element type Add/edit/delete a. Click on Add b. Select Structural LINK family elements c. Select 3D spar 8. d. Click on ok e. Click on close 4. Define Real constants Main menu: Preprocessor Real constants add/edit/delete a. Click on add 4
b. Click on OK c. Enter 0.5 for area d. Click on ok and close 5. Define material properties Main menu: preprocessorMaterial Properties Material models Structural Linear elastic isotropic a. Click on ok for material type 1 b. Enter 30e6 and 0.25 for PRXY c. Click on ok 6. Define nodes Main menu: Preprocessor Modelling_create Nodes In active CS a. Enter 1 for node number b. Enter (0,0,30) for X,Y,Z locations c. Enter 2 for node number d. Enter (0,0,-30) for X,Y,Z locations e. Enter 3 for node number f. Enter (0,-30,0) for X,Y,Z locations g. Enter 4 or node number h. Enter (40,0,0) for X,Y,Z locations 7. Define elements Main menu: Preprocessor modelling Create Elements Auto numbered Through nodes a. Pick the nodes 1 and 4. b. Click on apply c. Pick the nodes 2 and 4 d. Click on apply e. Pick the nodes 3 and 4 f. Click on apply
Solution 8. Apply displacements Main menu: Solution Define loads Apply Structural displacement On nodes a. Pick the node 1,2 and 3. b. Click on apply c. Select UX,UY,UZ d. Enter 0 against value e. Click on ok
5
9.Apply loads Main menu: Solution Loads Apply Structural force/moment on nodes a. Pick the node 4. b. Click on apply. c. Select FY . d. Enter -5000 for value. e. Click on ok. 10. Solve Min menu: Solution Solve Current LS a. click ok and close. Post Processing 11. Read results Main menu: General postprocRead results First set 12. Print deformed results Main menu: General postproc List results Nodal solution a. Select DOF solution. b. Select displacement sum c. Click on ok 13. Plot member forces Main menu: General postproc Element table Define table a. Add element table data b. Enter Labs as MEMBER c. Scroll data and choose sequence number d. Choose SMISC e. Enter 1 in the selection field after SMISC f. Click on ok g. Close 13. Print deformed shape Main menu:General postproc Plot results a. Click on deformed + undeformed b. Click ok 14. Exit Toolbar :quit a. Click on quit b. Click on ok.
6
Results The results of the matlab code agrees with the solution obtained by solving the problem in ANSYS.
The nodal displacements as obtained from the ansys and matlab code are tabled as follows. Node 1 2 3 4
UX 0 0 0 0.01736
UY 0 0 0 -0.06944
UZ 0 0 0 0
The member forces as obtained from the ansys and matlab code are tabled as follows Element no 1 2 3
Forces(lb) 4166.7 4166.7 -8333.33
Stresses(lb/in2) 8333.33 8333.33 -16666.66
Fig 1 The loading and boundary conditions in 3d truss
7
Fig 2 The deformed geometry of the 3D truss
Fig 3 The maximum and the minimum forces acting on the 3d truss
8
