added files for Undergraduate thesis's code

LYXSUPER · Nov 17, 2021 · 3d41a76 · 3d41a76
commit 3d41a76
Show file tree

Hide file tree

Showing 30 changed files with 1,588 additions and 0 deletions.
diff --git a/Assembling_Ksx_Ksy.m b/Assembling_Ksx_Ksy.m
@@ -0,0 +1,28 @@
+function [ Tsx, Tsy ] = Assembling_Ksx_Ksy
+global gcoord sdof
+Sx = find(gcoord(:,2)==max(gcoord(:,2))/2);  % at y = Hy/2 (along X-axes);
+Sy = find(gcoord(:,1)==max(gcoord(:,1))/2);  % at x = Lx/2 (along Y-axes);
+n = length(Sx) ;
+m = length(Sy) ;
+dofSx = zeros(1,2*n) ;
+dofSy = zeros(1,2*m) ;
+for i = 1:n
+    i1 = 2*i ;      i2 = i1-1 ;
+    dofSx(i1)=3*Sx(i);    dofSx(i2)=3*Sx(i)-2; 
+%     dw/dx                     w
+end   
+Tsx = sparse(2*n,sdof);
+for i=1:2*n
+    Tsx(i,dofSx(i))=1;
+end
+for j=1:m
+    j1 = 2*j ;      j2 = j1-1 ;
+    dofSy(j1)=3*Sy(j)-1;   dofSy(j2)=3*Sy(j)-2;
+%     -dw/dy                     w
+end
+Tsy = sparse(2*m,sdof);
+for j=1:2*m
+    Tsy(j,dofSy(j))=1;
+end
+end
+
diff --git a/Assembling_Stiffened.m b/Assembling_Stiffened.m
@@ -0,0 +1,11 @@
+function [Kp]=Assembling_Stiffened(option1,Emodule,Kp)
+if option1 == 1
+    [ Ksx ] = cal_Ksx(Emodule);
+    [ Tsx, Tsy ] = Assembling_Ksx_Ksy;
+    Kp=Kp+Tsx'*Ksx*Tsx;
+elseif option1 == 2
+    [ Ksx ] = cal_Ksx(Emodule);
+    [ Ksy ] = cal_Ksy(Emodule);
+    [ Tsx, Tsy ] = Assembling_Ksx_Ksy;
+    Kp = Kp + Tsx'*Ksx*Tsx + Tsy'*Ksy*Tsy ;
+end
diff --git a/Example_1.m b/Example_1.m
@@ -0,0 +1,9 @@
+Poisson=0.3;          % Poisson's ratio
+Emodule=1.7*10^7;  % Young elastic modulus  
+Load=1;      % Uniform load
+% Plate
+Lx=1;               % Length of the Plate along X-axes
+Hy=1;              % High of the Plate along Y-axes
+t=0.01;               % Thickness of plate
+% Beam // ox
+bx=0.01;    tx=0.1;     ISy= bx*tx^3/12;    %ex=tx/2;    Ax=bx*tx;   ISy= bx*tx^3/12;%+ex^2*Ax;
diff --git a/Example_2.m b/Example_2.m
@@ -0,0 +1,12 @@
+% Barick
+Poisson=0.3;          % Poisson's ratio
+Emodule=2.0684*10^5;  % Young elastic modulus  
+Load=6.89*10^-2;      % Uniform load
+% Plate
+Lx=762;               % Length of the Plate along X-axes
+Hy=1524;              % High of the Plate along Y-axes
+t=6.35;               % Thickness of plate
+% Beam // ox
+bx=12.7;    tx=127;     ISy= bx*tx^3/12;    %ex=tx/2;    Ax=bx*tx;   ISy= bx*tx^3/12;%+ex^2*Ax;
+% Beam // oy
+by=12.7;    ty=76.2;    ISx=by*ty^3/12;     %ey=ty/2;    Ay=by*ty;   ISx=by*ty^3/12;%+ey^2*Ay;
diff --git a/Grid_16_16_SP.fig b/Grid_16_16_SP.fig
diff --git a/Grid_32_32_SP.fig b/Grid_32_32_SP.fig
diff --git a/Grid_64_64.fig b/Grid_64_64.fig
diff --git a/Load_HCT_Sub_T1_SP.m b/Load_HCT_Sub_T1_SP.m
@@ -0,0 +1,30 @@
+function [fe]=Load_HCT_Sub_T1_SP(Load,x,y)
+x1=x(1);    x2=x(2);    x3=x(3);    y1=y(1);    y2=y(2);    y3=y(3);
+x0=(x1+x2+x3)/3;    y0=(y1+y2+y3)/3;
+% A12=2*area_Sub_triangle_1
+% o: The centroid of triangle element
+
+%                                          3-node (2)
+%                                            o
+%                                           / \
+
+%								L1=0      /		\   L3=0
+
+%										/ Sub-T_1 \
+
+%                         0-node  (3) o-------------o 2-node (1)
+%                                            L2=0 
+
+b30=y3-y0;                  b02=y0-y2;                  b23=y2-y3;
+c30=x3-x0;                  c02=x0-x2;                  c23=x2-x3;
+l30=sqrt(b30^2+c30^2);      l02=sqrt(b02^2+c02^2);      l23=sqrt(b23^2+c23^2);
+% mu30=(l23^2-l02^2)/l30^2;   mu02=(l30^2-l23^2)/l02^2;   
+mu23=(l02^2-l30^2)/l23^2;
+A12=b30*(-c02)-(-c30)*b02;      % A12=2*area_Sub_triangle_1
+A14=2*A12;
+
+fe1 =[ (5*b02*c23 - 5*b23*c02 - 2*b02*c30 + 2*b30*c02 + 7*b23*c30 - 7*b30*c23 - b02*c23*mu23 + b23*c02*mu23 - b23*c30*mu23 + b30*c23*mu23)/(80*(b02*c23 - b23*c02)), ((A12*b02*l23^2)/120 - (7*A12*b23*l23^2)/240 + (A12*b23*l23^2*mu23)/240)/(l23^2*(b02*c23 - b23*c02)) + (A14*c23)/(120*l23^2), (A14*b23)/(120*l23^2) - ((A12*c02*l23^2)/120 - (7*A12*c23*l23^2)/240 + (A12*c23*l23^2*mu23)/240)/(l23^2*(b02*c23 - b23*c02)), (7*b02*c23 - 7*b23*c02 - 2*b02*c30 + 2*b30*c02 + 5*b23*c30 - 5*b30*c23 + b02*c23*mu23 - b23*c02*mu23 + b23*c30*mu23 - b30*c23*mu23)/(80*(b23*c30 - b30*c23)), ((7*A12*b23*l23^2)/240 - (A12*b30*l23^2)/120 + (A12*b23*l23^2*mu23)/240)/(l23^2*(b23*c30 - b30*c23)) + (A14*c23)/(120*l23^2), (A14*b23)/(120*l23^2) - ((7*A12*c23*l23^2)/240 - (A12*c30*l23^2)/120 + (A12*c23*l23^2*mu23)/240)/(l23^2*(b23*c30 - b30*c23)), -(b02*c23 - b23*c02 - b02*c30 + b30*c02 + b23*c30 - b30*c23)/(20*(b02*c30 - b30*c02)), (A12*(b02 - b30))/(60*(b02*c30 - b30*c02)), -(A12*(c02 - c30))/(60*(b02*c30 - b30*c02))];
+
+fe=Load*A12*fe1';
+end
+
diff --git a/Load_HCT_Sub_T2_SP.m b/Load_HCT_Sub_T2_SP.m
@@ -0,0 +1,31 @@
+function [fe]=Load_HCT_Sub_T2_SP(Load,x,y)
+x1=x(1);    x2=x(2);    x3=x(3);    y1=y(1);    y2=y(2);    y3=y(3);
+x0=(x1+x2+x3)/3;    y0=(y1+y2+y3)/3;
+% A22=2*area_Sub_triangle_2
+% o: The centroid of triangle element
+
+%                                          3-node (1)
+%                                            o
+%                                           / \
+
+%								L3=0      /		\   L2=0
+
+%										/ Sub-T_2 \
+
+%                         1-node  (2) o-------------o 0-node (3)
+
+% x0=(x1+x2+x3)/3;   y0=(y1+y2+y3)/3;
+b10= y1-y0;                 b03=y0-y3;                  b31=y3-y1;
+c10= x1-x0;                 c03=x0-x3;                  c31=x3-x1;
+l10=sqrt(b10^2+c10^2);      l03=sqrt(b03^2+c03^2);      l31=sqrt(b31^2+c31^2);
+% mu10=(l31^2-l03^2)/l10^2;   mu03=(l10^2-l31^2)/l03^2;   
+mu31=(l03^2-l10^2)/l31^2;
+A22=b10*(-c03)-(-c10)*b03;      % A22=2*area_Sub_triangle_2
+A24=2*A22;
+
+
+fe2 =[ -(2*b03*c10 - 2*b10*c03 - 5*b03*c31 + 5*b31*c03 + 7*b10*c31 - 7*b31*c10 + b03*c31*mu31 - b31*c03*mu31 - b10*c31*mu31 + b31*c10*mu31)/(80*(b03*c31 - b31*c03)), ((A22*b03*l31^2)/120 - (7*A22*b31*l31^2)/240 + (A22*b31*l31^2*mu31)/240)/(l31^2*(b03*c31 - b31*c03)) + (A24*c31)/(120*l31^2), (A24*b31)/(120*l31^2) - ((A22*c03*l31^2)/120 - (7*A22*c31*l31^2)/240 + (A22*c31*l31^2*mu31)/240)/(l31^2*(b03*c31 - b31*c03)), (2*b03*c10 - 2*b10*c03 - 7*b03*c31 + 7*b31*c03 + 5*b10*c31 - 5*b31*c10 - b03*c31*mu31 + b31*c03*mu31 + b10*c31*mu31 - b31*c10*mu31)/(80*(b10*c31 - b31*c10)), (A24*c31)/(120*l31^2) - ((7*A22*b31*l31^2)/240 - (A22*b10*l31^2)/120 + (A22*b31*l31^2*mu31)/240)/(l31^2*(b10*c31 - b31*c10)), ((7*A22*c31*l31^2)/240 - (A22*c10*l31^2)/120 + (A22*c31*l31^2*mu31)/240)/(l31^2*(b10*c31 - b31*c10)) + (A24*b31)/(120*l31^2), (b03*c10 - b10*c03 - b03*c31 + b31*c03 + b10*c31 - b31*c10)/(20*(b03*c10 - b10*c03)), (A22*(b03 - b10))/(60*(b03*c10 - b10*c03)), -(A22*(c03 - c10))/(60*(b03*c10 - b10*c03))];
+
+fe=Load*A22*fe2';
+end
+
diff --git a/Load_HCT_Sub_T3_SP.m b/Load_HCT_Sub_T3_SP.m
@@ -0,0 +1,29 @@
+function [fe]=Load_HCT_Sub_T3_SP(Load,x,y)
+x1=x(1);    x2=x(2);    x3=x(3);    y1=y(1);    y2=y(2);    y3=y(3);
+x0=(x1+x2+x3)/3;    y0=(y1+y2+y3)/3;
+% o: The centroid of triangle element
+
+%                                          o-node (3)
+%                                            o
+%                                           / \
+
+%								L2=0      /		\   L1=0
+
+%										/ Sub-T_1 \
+
+%                         1-node  (1) o-------------o 3-node (2)
+%                                            L3=0 
+% x0=(x1+x2+x3)/3;   y0=(y1+y2+y3)/3;
+b20= y2-y0;                 b01=y0-y1;                  b12=y1-y2;
+c20= x2-x0;                 c01=x0-x1;                  c12=x1-x2;
+l20=sqrt(b20^2+c20^2);      l01=sqrt(b01^2+c01^2);      l12=sqrt(b12^2+c12^2);
+% mu20=(l12^2-l01^2)/l20^2;   mu01=(l20^2-l12^2)/l01^2;   
+mu12=(l01^2-l20^2)/l12^2;
+A32=b20*(-c01)-(-c20)*b01;      % A22=2*area_Sub_triangle_3
+A34=2*A32;
+
+fe3 =[ (5*b01*c12 - 5*b12*c01 - 2*b01*c20 + 2*b20*c01 + 7*b12*c20 - 7*b20*c12 - b01*c12*mu12 + b12*c01*mu12 - b12*c20*mu12 + b20*c12*mu12)/(80*(b01*c12 - b12*c01)), ((A32*b01*l12^2)/120 - (7*A32*b12*l12^2)/240 + (A32*b12*l12^2*mu12)/240)/(l12^2*(b01*c12 - b12*c01)) + (A34*c12)/(120*l12^2), (A34*b12)/(120*l12^2) - ((A32*c01*l12^2)/120 - (7*A32*c12*l12^2)/240 + (A32*c12*l12^2*mu12)/240)/(l12^2*(b01*c12 - b12*c01)), (7*b01*c12 - 7*b12*c01 - 2*b01*c20 + 2*b20*c01 + 5*b12*c20 - 5*b20*c12 + b01*c12*mu12 - b12*c01*mu12 + b12*c20*mu12 - b20*c12*mu12)/(80*(b12*c20 - b20*c12)), ((7*A32*b12*l12^2)/240 - (A32*b20*l12^2)/120 + (A32*b12*l12^2*mu12)/240)/(l12^2*(b12*c20 - b20*c12)) + (A34*c12)/(120*l12^2), (A34*b12)/(120*l12^2) - ((7*A32*c12*l12^2)/240 - (A32*c20*l12^2)/120 + (A32*c12*l12^2*mu12)/240)/(l12^2*(b12*c20 - b20*c12)), -(b01*c12 - b12*c01 - b01*c20 + b20*c01 + b12*c20 - b20*c12)/(20*(b01*c20 - b20*c01)), (A32*(b01 - b20))/(60*(b01*c20 - b20*c01)), -(A32*(c01 - c20))/(60*(b01*c20 - b20*c01))];
+
+fe=Load*A32*fe3';
+end
+
diff --git a/Main_Stiffened_Plate_HCT.m b/Main_Stiffened_Plate_HCT.m
@@ -0,0 +1,118 @@
+% Programming for Kirchhoff's plate bending using conforming HCT element
+% Static Analysis of  plate
+% Problem : To find the maximum bedning of plate when uniform transverse
+% pressure is applied. 
+% Boundary conditions is used, clamped
+%--------------------------------------------------------------------------
+% Code is written by : Tran Van Nha                                          |
+%                   Math & Compute Science                                |
+%                                                                         |
+% E-mail : [email protected]                                             |
+%--------------------------------------------------------------------------
+%|________________________________________________________________________|
+%|     Edge  :3Node [w,Qx,Qy]            Model:4Node-12Dof                |
+%|     Center:1Node [w]                                                   |
+%|                                                                        |
+%|________________________________________________________________________|
+%--------------------------------------------------------------------------
+% clc
+clear all
+tic
+format long
+global sdof nel nnel nnode Lx Hy t nx ny Poisson D Dmat
+global gcoord ele_nods H ndof edof ISy ISx
+%----------------------------------------------%
+%             1. PREPROCESSOR PHASE            %
+%----------------------------------------------%
+%------------------------------------------%
+%  1.1 Input the initial data of example   % 
+%------------------------------------------%
+% nx=input('Number of element in x-axis (8,12,16,20,24,...) = ');    
+% ny=input('Number of element in y-axis (8,12,16,20,24,...) = ');
+nx=6;                 % Number of element in x-axis
+ny=6;                 % Number of element in y-axis
+% % % % % % % % % % % % % % % % % % % % % % 
+fprintf('Choose example \n')
+fprintf('1. Single beam_ bxt \n')
+fprintf('2. Cross Stiffened_Barik \n')
+option1=input('Option (1 or 2) = ');
+
+if option1 == 1
+    Example_1;
+%     type_boundary
+elseif option1 == 2
+    Example_2;
+end
+
+%----------------------------------------------%
+%  1.2 Compute necessary data from input data  % 
+%----------------------------------------------%
+nel=2*nx*ny;          % Total number of elements in system
+nnel=3;               % Number of nodes per element
+ndof=3;               % Number of dofs per node
+edof=nnel*ndof;       % Degrees of freedom per element
+nnode=(nx+1)*(ny+1);  % Total number of nodes in system
+sdof=nnode*3;         % Total degrees of freedom in system
+D1=(Emodule*t^3)/(12*(1-Poisson^2));
+Dmat=[1 Poisson 0;Poisson 1 0;0 0 (1-Poisson)/2];  % Matrix of material constants
+H = D1*Dmat;                    % The Hooke matrix of Krichhoff's plate bending
+
+[gcoord,ele_nods,bcdof,bcval]=get_initdata_HCT_SP;
+%   gcoord     = matrix containing coordinate values of each node
+%   ele_nods   = matrix containing nodes of each element 
+%   bcdof      = vector containing dofs associated with boundary conditions
+%   bcval      = vector containing boundary condition values associated with dofs in bcdof
+
+Showmesh_Triangle(gcoord,ele_nods)
+
+%-----------------------------------------------%
+%             2. SOLUTION PHASE                 %
+%-----------------------------------------------%
+%-----------------------------------------------%
+%      2.1 Compute the stiffness matrix         % 
+%-----------------------------------------------%
+[ Kp,Fp ] = cal_Kp_Fp(Load);
+[Kp]=Assembling_Stiffened(option1,Emodule,Kp);
+% [ Ksx ] = cal_Ksx(Emodule);
+% % [ Ksy ] = cal_Ksy(Emodule);
+% % Assembling Ksx
+% [ Tsx, Tsy ] = Assembling_Ksx_Ksy;
+% % Kp = Kp + Tsx'*Ksx*Tsx + Tsy'*Ksy*Tsy ;
+% Kp=Kp+Tsx'*Ksx*Tsx;
+%-----------------------------------------------%
+%       2.2 Apply the boundary condition        %
+%-----------------------------------------------%
+[Kp Fp]=apply_bcdof_SP(Kp,Fp,bcdof,bcval);
+%-----------------------------------------------------------------------%
+% 2.3 Solve the system of equations to obtain nodal displacement vector %
+%-----------------------------------------------------------------------%
+% %-----------------------------------------------------------------------%
+% % 2.3 Solve the system of equations to obtain nodal displacement vector %
+% %-----------------------------------------------------------------------%
+% %-----------------------------------------------%
+% %              3. POSTPROCESSOR PHASE           %
+% %-----------------------------------------------%
+U=Kp\Fp;
+m=1:3:sdof;
+% Plate's deflection w:
+w=U(m);
+clear m;
+fh = figure ;
+set(fh,'name','Preprocessing for FEA','numbertitle','off','color','w') ;
+x=0:Lx/nx:Lx;
+y=0:Hy/ny:Hy;
+for j=1:ny+1
+    for i=1:nx+1
+        node=(j-1)*(nx+1)+i;
+        z(i,j)=U(node*3-2);
+    end
+end
+surf(x,y,z);
+title('Finite Element Mesh of Plate') ;
+% wmax=max(w)
+c1 = find(gcoord(:,2)==max(gcoord(:,2))/2);  % at y = Hy/2 (along Y-axes);
+c2 = find(gcoord(:,1)==max(gcoord(:,1))/2);  % at x = Lx/2 (along X-axes);
+cc=intersect(c1,c2);% taam
+wc_fem=U(cc*3-2)
+wcss=wc_fem*100*D1/(Load*(Lx)^4)
+toc