Essay: Finite element formulation

The finite element analysis includes the discretization of a system in its sub domain or
finite elements to solve a particular problem so as to achieve a greater accuracy and fast
analysis processes. To solve for the structural problems and various boundary value
problems there is a need of the finite element analysis. The finite element method can be
modeled for simple elements as well as for the complex shapes.
FEM can be used to solve the problems related to the heat conduction, fluid dynamics,
seepage flow, electric and magnetic field etc. As this method having a vast scope in
almost every engineering field, it prompted mathematicians to use this technique for the
solution of complicated boundary values.
 Consider the thermal equilibrium of an element of a heated 1-D body which is heated
from a side.
i) The rate at which heat enters from the left face say,
T
qx kA
x
ï,¶
 
ï,¶
…(57)
where, k is the thermal conductivity,
A is the cross sectional area through which heat flows (which is measured
perpendicular to the direction of heat flow)
(âˆ,T/âˆ,x) is the rate of change of temperature(T) w.r.t axial direction.
ii) The rate at which heat leaves the right face (by retaining the two terms of Taylor’s
series expansion)
. x
qx T T
q dx qx dx kA kA dx
x x x x
ï,¶  ï,¶ ï,¶  ï,¶  
         ï,¶  ï,¶ ï,¶  ï,¶  
…(58)
Now, energy balance for element for small time (dt)
Heat inflow in Heat generated by internal Heat outflow in Change in internal
Time (dt) + sources in time (dt) = time(dt) + energy during time(dt)
i.e ( . ) ( . ) x x dx
T
q dt qA dX dt q dt C dx dt
t
 
ï,¶
  
ï,¶
…(59)
where, q́ is the rate of heat generation per unit volume(by the heat source)
C is the specific heat , 𝜌 is the density
.
T
dt dT
t
ï,¶

ï,¶
(temperature change of the element in time dt)
T T
kA dx qA C
x x t

ï,¶  ï,¶  ï,¶
   
ï,¶  ï,¶  ï,¶
…(60)
From equation 60th we can see that the following cases could be formed such as:
 If the heat source q́ = 0, we get Fourier Equation,
T T
kA dx C
x x tIf the system is in steady state , we get Poission Equation,
T
kA dx qA
x x
ï,¶  ï,¶ 
   
ï,¶  ï,¶ 
…(62)
 If the heat source is zero and the system is in the steady state, we get the Laplace
Equation,
T
kA dx
x x
ï,¶  ï,¶ 
  
ï,¶  ï,¶ 
…(63)
 If the thermal conductivity and the area of cross section are constant. The equation
reduces to
T
x x
ï,¶  ï,¶ 
  
ï,¶  ï,¶ 
…(64)
The following types of equations such as Fourier equation, poissions equation and
Laplace equation can be found using the FEM on a 1D body considering the thermal
equilibrium. And similarly the FEM analysis could be done on various structures for
which the discretization of the structre could be done in a systematic and efficient way to
solve the boundary value problems.
There are 3 major categories of boundary value problems:
I. Equilibrium Problems —
i) These can be modeled to find the steady state displacement or stress
distribution in case of solid mechanics problem.
ii) It can be used to find temperature or heat flux distribution in case of heat
transfer models.
In fluid mechanics FEM can be used to find the pressure or velocity distribution.
II. Eigen Value Problem —
i) These are the extensions to the equilibrium problems in which some critical values of some parameters such as natural frequency various mode shapes and buckling loads are to be determined in addition to the corresponding steady state configuration.
4.2 Advantages of FEA
ï,· The physical properties, which are intractable and complex for any closed bound solution, can be analyzed by this method.
ï,· It can be used for any geometry (may be regular or irregular).
ï,· It can take care of any boundary conditions.
ï,· Material anisotropy and non-homogeneity can be catered without much difficulty.
ï,· It can take care of any type of loading conditions.
ï,· In this method approximations are confined to small sub domains.
ï,· In this method, the admissible functions are valid over the simple domain and have nothing to do with boundary, however simple or complex it may be.
ï,· Enable to computer programming.
4.3 Disadvantages of FEA
ï,· Computational time involved in the solution of the problem is high in case of complex problems.
For fluid dynamics problems some other methods of analysis may prove efficient than the FEM.
4.4 Limitations of FEA
ï,· Proper engineering judgment is to be exercised to interpret results.
ï,· It requires large computer memory and computational time to obtain intend results.
ï,· There are certain categories of problems where other methods are more effective, e.g., fluid problems having boundaries at infinity are better treated by the boundary element method.
ï,· For some problems, there may be a considerable amount of input data. It is always desirable to make a visual check of the input data.
ï,· In the FEM, many problems lead to round-off errors. Computer works with a limited number of digits and solving the problem with restricted number of digits may not yield the desired degree of accuracy or it may give total erroneous results in some cases. For many problems the increase in the number of digits for the purpose of calculation improves the accuracy.
4.5 Errors and Accuracy in FEA
Every physical problem is formulated by simplifying certain assumptions. Solution to the problem, classical or numerical, is to be viewed within the constraints imposed by these simplifications. The material may be assumed to be homogeneous and isotropic, its behaviour may be considered as linearly elastic, the prediction of the exact load in any
type of structure is next to impossible. As such the true behaviour of the structure is to be viewed with in these constraints and obvious errors creep in engineering calculations.
ï,· The results will be erroneous if any mistake occurs in the input data. As such, preparation of the input data should be made with great care.
ï,· When a continuum is discretised, an infinite degrees of freedom system is converted into a model having finite number of degrees of freedom. Thus the actual continuum is represented by a set of approximations.
ï,· The accuracy depends to a great extent on the mesh grading of the continuum. In regions of high precision, higher mesh grading is needed whereas in the other regions, the mesh chosen may be coarser. As the element size decreases, the discretisation error reduces.
ï,· Improper selection of shape of the element will lead to a considerable error in the solution. Triangle elements in the shape of an equilateral or rectangular element in the shape of a square will always perform better than those having unequal lengths of the sides. For very long shapes, the attainment of convergence is extremely slow.
ï,· In the finite element analysis, the boundary conditions are imposed at the nodes of the element whereas in an actual continuum, they are defined at the boundaries. Between the nodes, the actual boundary conditions will depend on the shape functions of the element forming the boundary.
ï,· Simplification of the boundary is another source of error. If the mesh is refined, then the error involved in the discretised boundary may be reduced.
MATLAB is a high-level language specially designed for dealing with matrices. This makes it particularly suited for programming the finite element method. In addition, MATLAB will allow the reader to focus on the finite element method by alleviating the programming burden.
ï,· Few Words on Writing MATLAB Programs
The Matlab programming language is useful in illustrating how to program the finite element method due to the fact it allows one to very quickly code numerical methods and has a vast predefined mathematical library. This is also due to the fact that matrix (sparse and dense), vector and many linear algebra tools are already defined and the developer can focus entirely on the implementation of the algorithm not defining these data structures.
The extensive mathematics and graphics functions further free the developer from the drudgery of developing these functions themselves or finding equivalent pre-existing libraries. A simple two dimensional finite element program in MATLAB need only be
few hundred lines of code whereas in FORTRAN or C++one might need a few thousand.
Although the Matlab programming language is very complete with respect to its mathematical functions there are a few finite element specific tasks that are helpful to develop as separate functions.
As usual there is a trade off to this ease of development. Since Matlab is an interpretive language each line of code is interpreted by the Matlab command line interpreter and executed sequentially at run time, the run times can be much greater than that of compiled programming languages like FORTRAN or C++.
It should be noted that the built-in Matlab functions are already compiled and are extremely efficient and should be used as much as possible. Keeping this slow down due to the interpretive nature of Matlab in mind, one programming construct that should be avoided at all costs is the for loop, especially nested for loops since these can make a MATLAB programs run time orders of magnitude longer than may be needed. Often for loops can be eliminated using Matlab’s vectorized addressing.
ï,· Algorithm of the solution
Before going to the result, it is important to see the general solution procedures to solve finite element problems using computer program. These procedures can be broadly classified into three sections:
1. Pre-processor – In this part of the program essential information is collected from the user about the geometry and type of the problem. Some of the major input parameters are listed below
ï,· Number of total nodes in the system
ï,· Number of total elements in the system
Coordinate values of every node in terms of the global coordinate system
ï,· Boundary conditions
ï,· Material property(Young’s Modulus and Poisson’s Ratio)
ï,· Specifying applied force
ï,· Size of the element i.e beam (thickness, length and width)
2. Processor — After taking user input all of the calculation will be conducted in this section of the program. Here major tasks of processor will be listed in this specific problem:
ï,· Calculating stiffness matrices and force vectors for every element
ï,· Assembling stiffness matrices into the system matrix and assembling force
vector to a system vector
ï,· Applying constraints (boundary conditions) to the system matrix and vector
ï,· Calculating deflection of the beam.
3. Postprocessor — Once the required calculation is performed in previous section one can display and plot the result in this part of the program.
ï,· Displaying the result, in this specific problem a deflection for each node will be displayed.
The eigen value problem for a beam element has been solved using FEM program written in Matlab and the results for the same has been compared to that of exact and the ABAQUS solution. The matlab files are used to compute both the element stiffness and mass matrices, so that they can be used to compute the natural frequency of a beam element. The element stiffness and mass matrices are assembled into system stiffness and system mass matrices. The matlab file feaplycs.m modifies the eigen value matrix
equation with the given constraints. This is used because instead of re-dimensioning the matrix size because of constraints, the original matrix size is conserved. The modified eigen value problem contain fictious zero eigen value in the same number of constraints. As a result the user can exclude these zero eigen values from the computer solution. In our problem we have seen that the solutions agree for the lower frequency. But the discrepencies become large for the higher modes of natural frequency. As a result of this we come to know that we can use more refined mesh to get the more accurate results for the higher modes to represent the corresponding mode shapes accurately. The refining of the mesh has been done both in our program that is the Beam Analyser (BA) as well as in the ABAQUS. And the results of both applications are compared with that of the Exact solutions.
The experimental procedure for cantilever beam are also modelled and the analysis has been performed using the impact hammer and the peak natural frequency for the element is found using the accelerometer and the FFT analyser.
ï,· Finite Element Method for Beam
 Problems
Problem 1: Simply supported beam under concentrated load at the centre.
Problem 2: Cantilever beam under free vibration using finite elements.
Problem 3: Experimental setup of a cantilever beam element using impact hammer and calculation of peak natural frequency.
Example Problem 5.2.1
Geometry
a) Cross-section of the Beam = (1 in * 1 in)
b) L/T ratio = 20
c) Length of Beam element = 20 in
Material Properties
d) Modulas of Elasticity (E) = (10 x 106 )psi
e) ν = 0.25
Boundary Condition
f) Simply Supported Beam with concentrated load of 100 lbs at center
 Example Problem 5.2.2
Geometry
a) Cross-section of the Beam = (0.02 m x 0.02 m)
b) L/T ratio = 50
c) Length of the Beam Element = 1 m
Material Properties
d) Modulas of Elastisity (E) = (100 x 109 )N/m2 and
e) Shear Modulas (G) = (40 x 109 )N/m2
f) ν = 0.25
Boundary Condition
g) Cantilever Beam
a) Cross-section of the Beam = (0.043 m x 0.005 m)
= (0.054 m x 0.005 m)
= (0.064 m x 0.005 m)
b) Length of the Beam Element (L1) = 0.460 m
(L2) = 0.550 m
(L3) = 0.660 m
Material Properties
c) Modulas of Elastisity (E) = (210 x 109 )N/m2 and
d) Shear Modulas (G) = (80 x 109 )N/m2
e) ν = 0.3
Boundary Condition
f) Cantilever Beam
5.3 Selecting the number of nodes per element
Simply if use the more number of nodes for a particular element the accuracy of the problem solution can be achieved but it will result in the greater computational efforts. So for this purpose a balance need to be maintained depending on the accuracy required and the time available for the solution. For a solution to be reached accurately the mesh size thereby can increased depending upon the problem size and the results can be computed close to the exact ones.
Based on the formulation described in the preceding sections, a finite element program called Beam Analyser (BA) is developed in MATLAB that can be used for linear static analysis and prediction of the behaviour of the beam element during the free vibration. The static analysis includes the slope and deflections caused in the beam and the free vibration includes the calculation of frequency and mode shapes.
For the purpose of validation of the Beam Analyser (BA) program and establishing its applicability, the above mentioned BA program is used to solve number of problems as mentioned above. The results of above problems obtained by finite element program are compared with the exact solutions and ABAQUS software. The beam element is divided into 12 number of elements or meshes in the beam vibration problem and the results are obtained for the same. The experimental results are obtained and the peak natural frequencies for the case of a cantilever beam is computed for varying length and also by varying the width of the beam and keeping the depth of the beam as constant.
Table 6.1: Different length for simply supported beam under concentrated load
at the centre:
Concentrated / Simply Supported
Slope(rad)
Length (in)
B.A
(Beam Analyser)
EXACT
(BA-EXACT)/ (BA)x100 % Diff
ABAQUS
(BA- ABAQUS)/ (BA)x100 % Diff
-0.0030
-0.0030
0.0
-0.0030
0.0
2
-0.0029
-0.0029
0.0
-0.0029
0.0
4
-0.0025
-0.0025
0.0
-0.0025
0.0
6
-0.0019
-0.0019
0.0
-0.0019
0.0
8
-0.0011
-0.0011
0.0
-0.0011
0.0
10
0.0000
0.0000
0.0
0.0000
0.0
12
0.0011
0.0011
0.0
0.0011
0.0
14
0.0019
0.0019
0.0
0.0019
0.0
16
0.0025
0.0025
0.0
0.0025
0.0
18
0.0029
0.0029
0.0
0.0029
0.0
20
-0.0030
-0.0030
0.0
-0.0030
0.0
Figure 6.1: Variation in slope through the length of a simply supported beam with concentrated load at the centre.
Table 6.2: Different Length for Simply supported beam under concentrated load
at the centre:
Concentrated /Simply supported
Deflection Z (m)
Length (in)
B.A
(Beam Analyser)
EXACT
(BA-EXACT)/ (BA)x100 % Diff
ABAQUS
(BA- ABAQUS)/ (BA)x100 % Diff
0.0000
0.0000
0.00
0.0000
0.0000
2
0.0059
0.0059
0.00
0.0059
0.0000
4
0.0114
0.0114
0.00
0.0114
0.0000
6
0.0158
0.0158
0.00
0.0159
-0.0063
8
0.0189
0.0189
0.00
0.0190
-0.0053
10
0.0200
0.0200
0.00
0.0201
-0.0050
12
0.0189
0.0189
0.00
0.0190
-0.0053
14
0.0158
0.0158
0.00
0.0159
-0.0063
16
0.0114
0.0114
0.00
0.0114
0.0000
18
0.0059
0.0059
0.00
0.0059
0.0000
20
0.0000
0.0000
0.00
0.0000
0.0000
Figure 6.2: Variation in deflection through the length of a simply supported beam with concentrated load at the centre.
0.0000
0.0050
0.0100
0.0150
0.0200
0.0250
26.2) Free Vibration Analysis (Eigen Value)
Table 6.3: Calculation of Different modes of free vibrations of a Cantilever Beam
Cantilever Beam
Mode 1 for 12 number of Elements
(b/d) Ratios
B.A
(Freq.)
EXACT
(Freq.)
(BA-Exact)/BA*100 %
ABAQUS
(Freq.)
(BA-ABAQUS)/BA*100 %
1
202.54
203.00
0.23
200.54
0.99
0.9
222.77
223.30
0.24
220.58
0.98
0.83
243.00
243.60
0.25
240.58
0.99
0.76
263.22
263.90
0.26
260.12
1.18
0.71
283.44
284.19
0.26
280.18
1.15
0.66
303.64
304.49
0.28
300.72
0.96
Figure 6.3: Variation in frequency for different (b/d) ratios of a cantilever beam.
Table 6.4: Calculation of Different modes of free vibrations of a Cantilever Beam
Cantilever Beam
Mode 2 for 12 number of Elements
(b/d) Ratios
B.A
(Freq.)
EXACT
(Freq.)
(BA-Exact)/BA*100 %
ABAQUS
(Freq.)
(BA-ABAQUS)/BA*100 %
1
1266.32
1272.17
0.46
1255.24
0.88
0.9
1392.09
1399.38
0.52
1380.31
0.85
0.83
1517.60
1526.60
0.59
1505.22
0.82
0.76
1642.86
1653.82
0.66
1630.01
0.78
0.71
1767.82
1781.03
0.74
1754.63
0.75
0.66
1892.48
1908.25
0.83
1879.02
0.71
Figure 6.4: Variation in frequency for different (b/d) ratios of a cantilever beam.
Table 6.5: Calculation of Frequency for different (b/d) ratios
Cantilever Beam
Mode 3 for 12 number of Elements
(b/d) Ratios
B.A
(Freq.)
EXACT
(Freq.)
(BA-Exact)/BA*100 %
ABAQUS
(Freq.)
(BA-ABAQUS)/BA*100 %
1
3553.10
3562.13
0.25
3528.85
0.68
0.9
3902.50
3918.34
0.40
3878.47
0.62
0.83
4250.25
4274.55
0.57
4227.17
0.54
0.76
4596.22
4630.77
0.75
4574.84
0.47
0.71
4940.28
4986.98
0.94
4921.48
0.38
0.66
5282.30
5343.19
1.14
5266.91
0.29
Figure 6.5: Variation in frequency for different (b/d) ratios of a cantilever beam
6.3) Experimental Analysis of Cantilever Beam
Table 6.6: Peak response of a cantilever beam using impact hammer and comparision of the result with B.A and ABAQUS
 Varying length of the element and constant width (0.043 m) and constant depth (0.005 m)
Length
(m) B.A Freq.(Hz) Experimental Freq. (Hz) (B.A-Experimental) / (B.A)*100% ABAQUS Freq.(Hz) (ABAQUS-Exp) /ABQ*100 (B.A-ABAQUS)/ B.A*100%
0.460 20.142 18.160 9 19.158 5 4
0.550 14.090 12.240 13 13.402 8 4
0.660 9.935 8.640 13 9.449 8 4
Figure 6.6: Variation in peak frequency for varying length of the beam element
Table 6.7: Peak response of a cantilever beam using impact hammer and comparision of the result with B.A and ABAQUS
 Varying length of the element and constant width (0.054 m) and constant depth (0.005 m)
Length
(m) B.A Freq.(Hz) Experimental Freq. (Hz) (B.A-Experimental) / (B.A)*100% ABAQUS Freq.(Hz) (ABAQUS-Exp) /ABQ*100 (B.A-ABAQUS)/ B.A*100%
0.460 20.142 17.840 11 19.158 6 4
0.550 14.090 12.080 14 13.402 9 4
0.660 9.785 8.320 14 9.307 10 4
Figure 6.7: Variation in peak frequency for varying length of the beam element
Table 6.8: Peak response of a cantilever beam using impact hammer and comparision of the result with B.A and ABAQUS
 Varying length of the element and constant width (0.064 m) and constant depth (0.005 m)
Length
(m) B.A Freq.(Hz) Experimental Freq. (Hz) (B.A-Experimental) / (B.A)*100% ABAQUS Freq.(Hz) (ABAQUS-Exp) /ABQ*100 (B.A-ABAQUS)/ B.A*100%
0.460 20.142 17.280 14 19.158 9 4
0.550 14.090 12.400 11 13.402 7 4
0.660 9.638 8.320 13 9.167 9 4
Figure 6.8: Variation in peak frequency for varying length of the beam element
It is evident that the Finite Element Method is a pretty robust means of numerically solving for vast majority of structural problems. The problem which we have considered is a beam element with varying support conditions and the same has been dealt for the calculation of slope and the deflection and the various parameters have been varied such as the number of elements used in the analysis. The free vibration studies are done by varying element size and and the higher modes of frequency are computed and are compared with that of the analytical solutions and the results obtained through ABAQUS. The observations made in this thesis can apply to even complex structures with minor changes like tailoring the right element and boundary condition.
Numerical results obtained by the Beam Analyser (B.A) are compared with that of the analytical solution and the results obtained through ABAQUS software.
 From Table 6.1 and 6.2 it may be concluded that the results shown are deviating not more than 1% than the exact and ABAQUS solution.
 The graphs have been plotted the slope and the deflection of a simply supported beam having a concentrated load at the centre is shown in Figure 6.1 and Figure 6.2.
 Form the Table 6.3 to Tble 6.5 the various mode frequencies are computed and the finite element processes have been used.
The BA program computes the various MODE frequencies for various (b/d) ratios of the beam element and shows a very good match with the exact solutions and the results have been validated through the ABAQUS software. By increasing the mesh size of the element and defining the number of degree of freedom required for the analysis one can easily compute for a more complex structure and for more number of elements having greater number of degree of freedom.
 The results of frequencies for various (b/d) ratios of the beam are shown in Figure 6.3 for MODE 1, Figure 6.4 for MODE 2 and in Figure 6.5 for MODE 3.
 The various mode shapes have been plotted from the Table 6.3 to Table 6.5 and they show a good convergence to the exact solutions and the ABAQUS results.
 The results for the experimental problem number 5.2.3 has been showed in the Table 6.6 for various length of the beam element and keeping the width as a constant factor for the particular 3 samples and it shows a good convergence with respect to the results obtained through BA and ABAQUS.
 In Table 6.7 and Table 6.8 the width of the beam element is varied (b = 0.054 & b = 0.064) keeping the length same as that in the above samples and the results are obtained and we see that the results also converges accordingly.
 As we increase the length and the width of the beam element we get that the peak natural frequencies gets decreased and as a result of these solutions various graphs have been plotted which can be seen if Figure 6.6 to Figure 6.8 which mainly shows the behaviour of the element if we vary the length and the width of the beam.
 It is also seen that the above solutions basically depends upon
Meshing
From Table 6.3 to Table 6.8 it can be seen that the result converges to the exact ones as we go on increasing the mesh size of the element. Meshing, is therefore an important aspect in running a numerical solution, and must always be treated with caution. A different mesh type can yield different answers depending on how the element code is set up.
Number of nodes per element
As we decide the size of the meshing which is to be used, the number of nodes per element governs the design.
1. Input Module
% Ex. 1
% Matlab program to solve a static beam deflection using hermitian beam elements
%
% variable descriptions
% k = element stiffness matrix
% kk = system stiffness matrix
% ff = system force vector
% index = a vector containing system dofs associated with each element
% bcdof = a vector containing dofs associated with boundary conditions
% bcval = a vector containing boundary condition values associated with the
% dofs in bcdof
%
nel = 5 ; % number of elements
nnel = 2 ; % number of nodes per element
ndof = 2 ; % number of dofs per node
nnode = (nnel-1)*nel + 1; % total number of nodes in system
sdof = nnode*ndof; % total system dofs
el = 10^7; % elastic modulas
xi = 1/12; % moment of inertia of crosssection
tleng = 10; % length of a half of the beam
leng = 10/nel;
area = 1 ; % crosssectional area of beam
rho = 1 ; % mass density(arbitary value for this problem
% because it is not used for the static problem)
ipt = 1; % option for mass matrix (arbitary values are not used here)
bcdof(1) = 1; % first dof (deflection at left end)is constrained
bcval(1) = 0; % whose described value is 0
bcdof(2) = 12; % 12th dof (slope at the symmetric end)is constrained
bcval(2) = 0; % whose described value is 0
ff = zeros(sdof,1); % initialization of system force vector
kk = zeros(sdof,sdof); % initialization of system matrix
index = zeros(nnel*ndof,1); % initialization of index vector
ff(11) = 50; % because a half of the load is applied due to symmetry
for iel = 1:nel % loop for total no. of elements
index = feeldof1(iel,nnel,ndof); % extract system dofs for each element
k = febeam1(el,xi,leng,area,rho,ipt); % compute element stiffness matrix
kk = feasmbl1(kk,k,index); % assembly into system matrix
end
[kk,ff] = feaplyc2(kk,ff,bcdof,bcval); % apply the boundary conditions
fsol = kkff; % solve the matrix equation
%
% analytical solution
%
e = 10^7 ; l = 20 ; xi = 1/12 ; P = 100 ;
for i = 1:nnode
x = (i-1)*2 ;
c = P/(48 * e * xi);
k = (i-1) * ndof + 1;
esol(k) = c*(3*l^2-4*x^2)*x;
esol(k+1) = c*(3*l^2-12*x^2);
end
%
% Print both exact and fem solutions
%
num = 1:1:sdof;
store = [num’ fsol esol’]
%
% Ex.2. MATLAB program to solve the natural frequencies of a beam using beam Elements with displacement degrees of freedom only.
%
% Variable Descriptions
% k = element stiffness matrix
% m = element mass matrix
% kk = system stiffness matrix
% mm = system mass matrix
% index = a vector containing system dofs associated with each element
% bcdof = a vector containing dofs associated with boundary conditions
%
nel = 12; % no of elements
nnel = 2 ; % no of nodes per element
ndof = 3 ; % no of dofs per node
nnode = (nnel-1)*nel+1 ; % total no. of nodes in the elemnt
sdof = nnode*ndof ; % total sysyem dofs
el = 100*10^9 ; % elastic modulas
sh = 40*10^9 ; % shear modulas
tleng = 1 ; % total length of the beam
leng = tleng/nel ; % uniform mesh (element size of beam)
heig = 0.03 % height (or thickness)of the beam
width = 0.02 % width of the beam
rho = 1000 ; % mass density of the beam
bcdof(1) = 1; % bottom inplane displacement at node 1 is constrained
bcdof(2) = 2; % top inplane displacement at node 1 is constrained
bcdof(3) = 3; % transverse disp. at node 1 is constrained
kk = zeros(sdof,sdof); % initialization of sys.
stiffness matrix
mm = zeros(sdof,sdof); % initialization of sys. mass matrix
index = zeros(nnel*ndof,1); % initialization of index vectors
for iel = 1:nel % loop for the total no. elements
index = feeldof1(iel,nnel,ndof); % extract system dofs for each element
[k,m] = febeam3(el,sh,leng,heig,width,rho);% compute element matrices
kk = feasmbl1(kk,k,index); % assembly of sys. stiffness matrix
mm = feasmbl1(mm,m,index); % assembly of sys. mass matrix
end
[kk,mm] = feaplycs(kk,mm,bcdof); % apply the boumdary conditions
fsol = eig (kk,mm); % solve the matrix equation and print
fsol = sqrt(fsol)
%
fsol = eig(kk,mm); % solve the eigen value problem
fsol = double(sqrt(fsol)) % print circular frequencies)
%
2. Function feeldof1
function [index] = feeldof1(iel,nnel,ndof)
%
% purpose
% compute system dofs associated with each element in one-D problem
%
% synopsis
% [index] = feeldof1(iel,nnel,ndof)
%
% variable description:
% index – system dof vector associated with element iel
% iel – element no. whose system dofs are to be determined
% nnel – no. of nodes per element
% ndof – no. of dofs per node
%
%
edof = nnel*ndof;
start = (iel-1)*(nnel-1)*ndof ;
%
for i = 1:edof
index(i) = start + i ;
end
%
3. Function febeam1
function [k,m] = febeam1(el,xi,leng,area,rho,ipt)
%
% purpose:
% [k,m] = febeam(el,xi,leng,area,rho,ipt)
%
% Variable description
% k – element stiffness matrix (size of 4*4)
% m – element mass matrix (size of 4*4)
% el – elastic modulas
% xi – second moment of inertia of cross section
% leng – element length
% area – area of beam cross section
% rho – mass density (mass per unit volume)
% ipt – 1:consistent mass matrix
% 2 : lumped mass matrix
% otherwise : diagonal mass matrix
%
%
% stiffness matrix
%
c = el*xi/(leng^3);
k = c*[12 6*leng -12 6*leng ;…
6*leng 4*leng^2 -6*leng 2*leng^2;…
-12 -6*leng 12 -6*leng;…
6*leng 2*leng^2 -6*leng 4*leng^2];
%
% consistent mass matrix
%
if ipt == 1
%
mm = rho*area*leng/420 ;
m = mm*[156 22*leng 54 -13*leng;…
22*leng 4*leng^2 13*leng -3*leng^2;…
54 13*leng 156 -22*leng;…
-13*leng -3*leng^2 -22*leng 4*leng^2];
%
% lumped mass matrix
%
elseif ipt == 2
%
m = zeros (4,4);
mass = rho*area*leng;
m = diag([mass/2 0 mass/2 o]);
%
% diagoal mass matrix
else
%
m = zeros (4,4);
mass = rho*area*leng;
m = mass*diag([1/2 leng^2/78 1/2 leng^2/78]);
%
end
%
4. Function febeam3
function [k,m] = febeam3(el,sh,leng,heig,width,rho)
%
% purpose:
% stiffness and mass matrices for beam element with displacement degree of freedom only
% nodal dof ub1 ut1 v1 ub2 ut2 v2
%
% synopsis
% [k,m] = febeam1(el,sh,leng,heig,rho,area,ipt)
%
% variable description:
% k – element stiffness matrix (size of 6*6)
% m – element mass matrix (size of 6*6)
% el – elastic modulas
% sh – shear modulas
% leng – element length
% heig – element thickness
% width – width of the beam element
% rho – mass density of the beam element(mass per unit volume)
% lumped mass matrix only
%
% stiffness matrix
%
a1 = (sh*leng*width)/(4*heig);
a2 = (sh*heig*width)/leng;
a3 = (el*heig*width)/(6*leng);
a4 = sh*width/2;
k = [ a1+2*a3 -a1+a3 a4 a1-2*a3 -a1-a3 -a4;…
-a1+a3 a1+2*a3 -a4 -a1-a3 a1-2*a3 a4;…
a4 -a4 a2 a4 -a4 -a2;…
a1-2*a3 -a1-a3 a4 a1+2*a3 -a1+a3 -a4;…
-a1-a3 a1-2*a3 -a4 -a1+a3 a1+2*a3 a4;…
-a4 a4 -a2 -a4 a4 a2];
%
% lumped mass matrix
%
m = zeros(6,6);
mass = rho*heig*width*leng/4;
m = mass*diag([1 1 2 1 1 2]);
%
5. Function feasmbl1
function [kk] = feasmbl1(kk, k, index)
edof = length(index);
for i=1:edof
ii = index(i);
for j = 1:edof
jj = index(j);
kk(ii,jj) = kk(ii,jj) + k(i,j);
end
end
6. Function feasmbl1
function [kk] = feasmbl1(mm, m, index)
edof = length(index);
for i=1:edof
ii = index(i);
for j = 1:edof
jj = index(j);
kk(ii,jj) = kk(ii,jj) + k(i,j);
end
end
7. Function feaplycs
function [kk,mm]= feaplycs(kk,mm,bcdof)
% Purpose:
% Apply constraints to eigenvalue matrix equation [kk]
% Synopsis:
% [kk,mm]= feaplycs(kk,mm,bcdof)
%
% Variable Description:
% kk – system stiffness matrix before applying constraints
% mm – system mass matrix before applying constraints
% bcdof – a vector containging constrained d.o.f
%
n = length(bcdof);
sdof = size(kk);
for i=1:n
c = bcdof(i);
for j=1:sdof
kk(c,j)=0;
kk(j,c)=0;
mm(c,j)=0;
mm(j,c)=0;
end
mm(c,c)=1;
end
%
8. Function feaplyc2
function [kk,ff] = feaplyc2(kk,ff,bcdof,bcval)
%
% purpose
% apply constraints to matrix equation [kk]x=ff
%
% synopsis
% [kk,ff] = feaplybc(kk,ff,bcdof,bcval)
%
% variable description
% kk – sysyem matrix before applying constraints
% ff – sysyem vector before applying constraints
% bcdof – a vector containing constraind dof
% bcval – a vector containing constrained value
%
% For example, there are constraints at dof = 2 and 10
% and their constrained values are 0.0 and 2.5, respectively.
% Then, bcdof(1)=2 and bcdof(2)=10;
% and bcval(2) = 2.5.
%
%
n = length(bcdof);
sdof = size(kk);
%
for i = 1:n
c = bcdof(i);
for j = 1:sdof
kk (c,j) = 0;
end
%
kk(c,c) = 1;
ff(c) = bcval(i);
end
%