-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathPlateBending.m
176 lines (159 loc) · 4.99 KB
/
PlateBending.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
%This program performs plate bending
% problems using C1 rectangular elements
%Functions will work on Octave, FreeMat
% and Matlab
%Created by Mohammad Tawfik
%mohammad.tawfik@gmail.com
%In assotiation with research papers
% published on ResearchGate.Net
%Author: Mohammad Tawfik
%Title: Fundamentals of Numerical Analysis
%DOI: 10.13140/RG.2.2.25680.81925
%Updated text link:
%https://www.researchgate.net/publication/321850359_Fundamentals_of_Numerical_Analysis_Book_Draft
%Title: Finite Element Analysis
%DOI: 10.13140/RG.2.2.32391.70560
%Updated text link:
%https://www.researchgate.net/publication/321850256_Finite_Element_Analysis_Book_Draft
%Title: Dynamics and Control of Flexible Structures
%DOI: 10.13140/RG.2.2.29036.26242
%Updated text link:
%https://www.researchgate.net/publication/321850001_Dynamics_and_Control_of_Flexible_Structures_Book_Draft
%Title: Panel Flutter
%DOI: 10.13140/RG.2.1.1537.6807
%Updated text link:
%https://www.researchgate.net/publication/275712979_Panel_Flutter
%More code about other topics in the text
% may be downloaded from:
% https://github.com/mohammadtawfik/FEM-Plates
%Clearing the memory
clear all
close all
clc
%Problem Data
%Aluminum material properties are used
Modulus = 71e9; %GPa - Modulus of elasticity
Nu = 0.3; % Poisson's ration
Rho = 2700; %Kg/m^3 - Density
Alpha = 22.5e-6;% Thermal expansion coefficient
%Geometric Data
LengthX = 1; % m
LengthY = 1; % m
Thickness= 0.001; % m
%Finite Element Problem Data
Nx = 10; %number of elements in the x-direction
Ny = 10; %number of elements in the x-direction
%Evaluating the basic quantities
Lx=LengthX/Nx; %element dimensions
Ly=LengthY/Ny;
Qq=Modulus/(1-Nu*Nu)* ...
[1,Nu,0;Nu,1,0;0,0,(1-Nu)/2];
Dd=Qq*Thickness*Thickness*Thickness/12;
TbInv=CalcTbInv(Lx,Ly);
[Kb,Mb]=CalcLinear(Dd,Lx,Ly);
Kb=TbInv'*Kb*TbInv;
Mb=TbInv'*Mb*TbInv*Rho*Thickness;
%***************************************
%Creating the Mesh (Noedes and Elements)
%***************************************
%Initializing the Nodes registry!
% Columns 1&2 contain the
% coordinates of the nodes (x,y)
% Columns 3 to 6 contain the numbers
% of the degrees of freedom at the node
Nodes=zeros((Nx+1)*(Ny+1),6);
%filling the Nodes registry
for jj=1:Ny+1
for ii=1:Nx+1
Pointer=(jj-1)*(Nx+1) + ii; %Node number
Nodes(Pointer,1)=(ii-1)*Lx;
Nodes(Pointer,2)=(jj-1)*Ly;
Nodes(Pointer,3)=Pointer*4-3; %w
Nodes(Pointer,4)=Pointer*4-2; %wx
Nodes(Pointer,5)=Pointer*4-1; %wy
Nodes(Pointer,6)=Pointer*4; %wxy
endfor
endfor
%Initializing the Elements registry
% Columns 1 to 4 contain the global node
% numbers associated with the element
Elements=zeros(Nx*Ny,4);
for jj=1:Ny
for ii=1:Nx
Pointer=(jj-1)*Nx + ii;
Elements(Pointer,1)=(jj-1)*(Nx+1)+ii;
Elements(Pointer,2)=(jj-1)*(Nx+1)+ii+1;
Elements(Pointer,3)=(jj)*(Nx+1)+ii+1;
Elements(Pointer,4)=(jj)*(Nx+1)+ii;
endfor
endfor
%
%Creating the Boundary Conditions vector
% This vector will contain the numbers of
% the fixed degrees of freedom;
% ie. the columns and rows that will be
% eleminated from the global matrix
BCs=[];
%Support conditions are user-defined
% 0=CCCC
% 1=SSSS
SupportConditions=0;
for jj=1:Ny+1
for ii=1:Nx+1
Pointer=(jj-1)*(Nx+1) + ii; %Node number
%If the node is on the sides ...
if or(or(jj==1,jj==Ny+1),or(ii==1,ii==Nx+1))
if SupportConditions==0
%The following line counts the BCs
% for a CCCC plate
BCs=[BCs,Nodes(Pointer,3:6)];
elseif SupportConditions==1
%The following lines counts for BCs
% for SSSS plate
%if the node is a corner node
if and(or(jj==1,jj==Ny+1),or(ii==1,ii==Nx+1))
BCs=[BCs,Nodes(Pointer,3:6)];
%if the node is on either horizontal sides
elseif or(jj==1,jj==Ny+1)
BCs=[BCs,Nodes(Pointer,3:4)];
%if the node is on either vertical sides
elseif or(ii==1,ii==Nx+1)
BCs=[BCs,Nodes(Pointer,3)];
BCs=[BCs,Nodes(Pointer,5)];
endif
endif
endif
endfor
endfor
%The complementary coundary conditions
% is a vector that contains the numbers
% of the degrees of freedom that are FREE
BCsC=1:4*(Nx+1)*(Ny+1);
%BCs
BCsC(BCs)=[];
%*****************************************
%The global Matrices
%*****************************************
%Initialization
KG=zeros(4*(Nx+1)*(Ny+1),4*(Nx+1)*(Ny+1));
MG=zeros(4*(Nx+1)*(Ny+1),4*(Nx+1)*(Ny+1));
%Looping for the elements
for ii=1:Nx*Ny
%Looping for the element nodes
for jj=1:4
DOFsj=Nodes(Elements(ii,jj),3:6);
%Looping AGAIN for the elements nodes;
for kk=1:4
DOFsk=Nodes(Elements(ii,kk),3:6);
KG(DOFsj,DOFsk)=KG(DOFsj,DOFsk)+ ...
Kb(4*jj-3:4*jj,4*kk-3:4*kk);
MG(DOFsj,DOFsk)=MG(DOFsj,DOFsk)+ ...
Mb(4*jj-3:4*jj,4*kk-3:4*kk);
endfor
endfor
endfor
%Applying the boundary conditions
KReduced=KG(BCsC,BCsC);
MReduced=MG(BCsC,BCsC);
vvv=sqrt(sort(eig(inv(MReduced)*KReduced)));
vvv(1:10)*sqrt(Rho*Thickness/Dd(1,1))