%% Deflection of a beam using a 4 node isoparametric element % Coded by(Matthias Schubert, 2007, PhD Seminar 'The Finite Element Method and the Analysis of Systems % with Uncertain Properties') clc; clear all; clf; %% Input %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% % Geometrical input l=1000; % Length [mm] h=60; % Height [mm] b_el=40; % Width [mm] Num_el_x=30; % Number of elements in x-direction Num_el_y=1; % Number of elements in y-direction Upper_nodes=1:(Num_el_y+1):(Num_el_x+1)*(Num_el_y+1)-Num_el_y; Lower_nodes=(Num_el_y+1):(Num_el_y+1):(Num_el_x+1)*(Num_el_y+1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% %Input forces % Force at node number F_Nodes_x=[];%[1+(Num_el_y+1)*Num_el_x]; F_Nodes_y=[1+(Num_el_y+1)*Num_el_x]; % Force Force_x=[]; %[N] Force_y=[-1000]; %[N] % Boundary conditions at nodes B_Node_x=[1 2 3]; %Node Number B_Node_y=[1 2 3]; %Node Number %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% % Material Properties nu=0.3; %Poisson ratio E=210000; %Stiffness [N/mm^2] % Aplification for plot amp=100; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Main Programm %% % Number of Elements N_El=Num_el_x*Num_el_y; %% % Number of nodes N_Nodes=(Num_el_x+1)*(Num_el_y+1); %% % Generating the Load Vector Load_glob=zeros(N_Nodes*2,1); for i=1:size(F_Nodes_y,2) Load_glob(F_Nodes_y(i)*2,:)=Force_y(i); end for i=1:size(F_Nodes_x,2) Load_glob(F_Nodes_x(i)*2-1,:)=Force_x(i); end %% % Generating boundary vector for i=1:size(B_Node_x,2) Boundary(i,1)=B_Node_x(i)*2-1; end for j=1:size(B_Node_y,2) Boundary(j+i,1)=B_Node_y(j)*2; end Boundary=sort(Boundary); %% % Generating nodal coordinates X_coordinates=zeros(Num_el_y+1,Num_el_x+1); Y_coordinates=zeros(Num_el_y+1,Num_el_x+1); for x_coor=1:(Num_el_y+1) X_coordinates(x_coor,:)=linspace(0,l,(Num_el_x+1)); end for y_coor=1:(Num_el_x+1) Y_coordinates(:,y_coor)=linspace(0,h,(Num_el_y+1)); end trans1=(1:1:Num_el_y+1)'; trans2=(size(trans1,1)+1)*ones(Num_el_y+1,1); trans=trans2-trans1; clear trans1 trans2; Y_coordinates=Y_coordinates(trans,:); %% % Assembling elements Index=[2 4 3 1]; ele_counter=0; Node_tmp_1=1; Node_tmp_2=Num_el_y+2; for m=1:Num_el_x for n=1:Num_el_y ele_counter=ele_counter+1; eval(sprintf('Element_%i=[Node_tmp_1 Node_tmp_1+1 Node_tmp_2 Node_tmp_2+1];',ele_counter)); eval(sprintf('Element_%i=Element_%i(Index);',ele_counter,ele_counter)); Node_tmp_1=Node_tmp_1+1; Node_tmp_2=Node_tmp_2+1; end Node_tmp_1=Node_tmp_1+1; Node_tmp_2=Node_tmp_2+1; end %% % Local element stiffness matrix (analytically derived) a_el=l/Num_el_x; h_el=h/Num_el_y; gamma=a_el/h_el; Psi1=(1+nu)*gamma; Psi2=(1-3*nu)*gamma; Psi3=2+(1-nu)*gamma^2; Psi4=2*gamma^2+(1-nu); Psi5=(1-nu)*gamma^2-4; Psi6=(1-nu)*gamma^2-1; Psi7=4*gamma^2-(1-nu); Psi8=gamma^2-(1-nu); k_4node =(E*b_el)/(24*gamma*(1-nu^2))*[4*Psi3 3*Psi1 2*Psi5 -3*Psi2 -2*Psi3 -3*Psi1 -4*Psi6 3*Psi2 ;... 3*Psi1 4*Psi4 3*Psi2 4*Psi8 -3*Psi1 -2*Psi4 -3*Psi2 -2*Psi7 ;... 2*Psi5 3*Psi2 4*Psi3 -3*Psi1 -4*Psi6 -3*Psi2 -2*Psi3 3*Psi1 ;... -3*Psi2 4*Psi8 -3*Psi1 4*Psi4 3*Psi2 -2*Psi7 3*Psi1 -2*Psi4 ;... -2*Psi3 -3*Psi1 -4*Psi6 3*Psi2 4*Psi3 3*Psi1 2*Psi5 -3*Psi2 ;... -3*Psi1 -2*Psi4 -3*Psi2 -2*Psi7 3*Psi1 4*Psi4 3*Psi2 4*Psi8 ;... -4*Psi6 -3*Psi2 -2*Psi3 3*Psi1 2*Psi5 3*Psi2 4*Psi3 -3*Psi1 ;... 3*Psi2 -2*Psi7 3*Psi1 -2*Psi4 -3*Psi2 4*Psi8 -3*Psi1 4*Psi4 ]; %% % Assemble global stiffness matrix K_global=sparse(N_Nodes*2,N_Nodes*2); for i=1:N_El for j=1:4 for k=1:4 eval(sprintf('K_global( Element_%i(1,j)*2-1, Element_%i(1,k)*2-1 )= K_global( Element_%i(1,j)*2-1, Element_%i(1,k)*2-1 ) + k_4node(j*2-1,k*2-1);',i,i,i,i)); eval(sprintf('K_global( Element_%i(1,j)*2-1, Element_%i(1,k)*2 )= K_global( Element_%i(1,j)*2-1, Element_%i(1,k)*2 ) + k_4node(j*2-1,k*2);',i,i,i,i)); eval(sprintf('K_global( Element_%i(1,j)*2, Element_%i(1,k)*2-1 )= K_global( Element_%i(1,j)*2, Element_%i(1,k)*2-1 ) + k_4node(j*2,k*2-1);',i,i,i,i)); eval(sprintf('K_global( Element_%i(1,j)*2, Element_%i(1,k)*2 )= K_global( Element_%i(1,j)*2, Element_%i(1,k)*2 ) + k_4node(j*2,k*2);',i,i,i,i)); end end end %% % Result K_global_tmp=K_global; K_global(:,Boundary)=[]; K_global(Boundary,:)=[]; Load_glob(Boundary,:)=[]; UG_result=K_global\Load_glob; %% % Creating deformation marix U=(zeros(N_Nodes*2,1)); Ug_init=(1:1:N_Nodes*2)'; Ug_init(Boundary,:)=[]; U(Ug_init,:)=UG_result; U_x_tmp=(U(1:2:end).*amp)'; U_y_tmp=(U(2:2:end).*amp)'; U_x=zeros(Num_el_y+1,Num_el_x+1); U_y=zeros(Num_el_y+1,Num_el_x+1); for t=1:(Num_el_y+1) eval(sprintf('U_x(%i,:)=U_x_tmp(%i:Num_el_y+1:end);',t,t)); eval(sprintf('U_y(%i,:)=U_y_tmp(%i:Num_el_y+1:end);',t,t)); end Deformation_x=X_coordinates+U_x; Deformation_y=Y_coordinates+U_y; %% Calculating stresses %Strain-displacement matrix B_mat=1/(2*a_el*h_el).*[ -h_el 0 h_el 0 h_el 0 -h_el 0 ;... 0 -a_el 0 -a_el 0 a_el 0 a_el ;... -a_el -h_el -a_el h_el a_el h_el a_el -h_el ]; %Stress-strain matrix E_mat=(E)/(1-nu^2).*[ 1 nu 0 ;... nu 1 0 ;... 0 0 0.5*(1-nu) ]; %Local displacement matrix for j=1:N_El for k=1:4 eval(sprintf('Element_%i_disp(k*2-1,1) =U(Element_%i(1,k)*2-1,1);',j,j)); eval(sprintf('Element_%i_disp(k*2,1) =U(Element_%i(1,k)*2,1);',j,j)); end eval(sprintf('Stress_mat_test(:,j)=B_mat*Element_%i_disp;',j)); eval(sprintf('Stress_mat(:,j)=E_mat*B_mat*Element_%i_disp;',j)); eval(sprintf('clear Element_%i_disp Element_%i',j,j)); end ele_counter=0; for k=1:Num_el_x for j=1:Num_el_y ele_counter=ele_counter+1; Sigma_xx(j,k)=Stress_mat(1,ele_counter); Sigma_yy(j,k)=Stress_mat(2,ele_counter); Sigma_xy(j,k)=Stress_mat(3,ele_counter); end end %% Plot results figure(1) AXIS equal GRID ON hold on TITLE('Deformed System','FontSize',18); %% % Plot undeformed system for coor=1:size(Y_coordinates,1) plot(X_coordinates,Y_coordinates(coor,:),'-r'); end for coor=1:size(X_coordinates,2) plot(X_coordinates(:,coor),Y_coordinates,'-r'); end %% % Plot deformed system plot(Deformation_x,Deformation_y,'-b'); for coor=1:Num_el_y+1 eval(sprintf('plot(Deformation_x(%i,:),Deformation_y(%i,:));',coor,coor)); end hold off %Plot stresses figure(2) subplot(3,1,1); imagesc(Sigma_xx); COLORBAR('EastOutside'); TITLE('Sigma xx [MPa]','FontSize',18); subplot(3,1,2); imagesc(Sigma_yy); COLORBAR('EastOutside'); TITLE('Sigma yy [MPa]','FontSize',18); subplot(3,1,3); imagesc(Sigma_xy); COLORBAR('EastOutside'); TITLE('Sigma xy [MPa]','FontSize',18); %% % Deleting Variables %clear Node_tmp* Psi* ele_counter clear i k t j coor U_y_tmp U_x_tmp Ug_init m n trans UG_result U_x U_y x_coor y_coor %% Stresses