Unformatted Attachment Preview
THEORETICAL BACKGROUND
In this project, the aim is to assume linear elasticity between the concrete and the
soil underneath the gravity dam. Fluids present will be taken as pressure loads and one
should also assume that the soil and concrete are perfectly bonded to each other.
Therefore, an analysis that’s meant to break down these data into simpler parts, so that
it can be reassembled with a supported system of equations to address the problem need
to conduct a Finite Element Analysis. It observes a Partial Differential Equations model
by nature, and an iterative process.
This method allows easier modeling and simulation of complex geometrical and
irregular shapes. It allows virtual testing in the sense that specifications can be adjusted
should physical properties of a certain project in a design process collide or do not match.
Boundary conditions can define which conditions need a response from FEM. Meshes
that are generated by the soft wares given for this activities consider different types of
elements. Although the problem does not really expect perfect shapes of the figures,
assumption of them being geometrically perfect shall be made to support computergenerated calculations done from Abaqus and as well as at the MATLAB.
Studying the theoretical background of the Brachistochrone Problem (“Shortest
Time”) would be by Galilei back in 1638 who discovered in his last work that a sphere
could travel faster from time A to B following that there is a curvature path followed, and
the distance traveled also has gravity influencing it. The problem is related with time and
delay, and it was only in 1696 when Bernoulli, the father of Variational calculus, was able
to solve this due to Jacob’s influence as Bernoulli shows his 15th version in Leipzig. Other
scientists and breakthroughs followed, such as from Leibniz in 1697 where he first drafted
the discrete variational method with elementwise triangular shape functions; Lagrange
where an arbitrary variation, delta y, is now the test function; Courant who introduced
triangular and rectangular “finite elements” for the 2DSt.Venant torsion problem of a
prismatic bar (Poisson equation); Isoparametric mapping for quadrilateral elements in
1968; model hierarchies in 1991 that separate errors of discretization and idealization,
which is essential for verification, validation, and uncertainty quantification. Until the
recent breakthroughs of availability of simulation for FEM which was initially introduced
by SimScale in 2013. Many people through time worked on improving this technique of
solving problems related to engineering and mathematical models not only because of
the method it offers but also one can view the techniques done via simulation.
Governing Equations of Gravity Dam
The DamReservoir interaction is represented by two coupled differential
equations of the second order. The equation of the structure and the reservoir can be
written in the form:
(1)
(2)
Where [M], [C] and [K] are mass, damping and stiffness matrices of the structure, and [G],
[C’] and [K’] are matrices representing mass, damping and stiffness of the reservoir, respectively.
[Q] is the coupling matrices and {f1} is the vector of body force and hydrostatic force. {f2} is the
component of the force due to acceleration at the boundaries of damreservoir and reservoir
foundation. {P} and {u} are the vector of pressure and displacement.
is the ground
acceleration and ρ is the density of the fluid. The dot represents the time derivative.
Plain Strain
Plain strain refers to the physical deformation of a body that is characterized by the
displacement of material in a direction that is parallel to the given plane. Gravity dam is a solid
plain strain structure. Its thickness is much greater than its other two dimensions, that is why it
has been analyzed as 2D plain strain structure. Gravity dams are subjected to various forces like
hydrostatic pressure, uplift pressure etc. due to which it causes stress concentration within its
body. Such stress concentration leads to localized failure zones in the structure. Though the
stress concentration is to be localized can leads to drastic damage to important structures like
dams. The dam structure failure is often analyzed using tools like Finite Element Method &
ABAQUS.
NUMERICAL (FEM) BACKGROUND:
For this section, the breakdown of some concepts and operations mostly used in
this report will be elaborated, such as the derivation of FEM weak forms, Galerkin's
method, and discretization by triangular elements, and post processing.
Divide and Conquer Approach
In the FEM, a continuous domain is discretized into simple geometric shapes called
elements. A finite element is a small piece of structure. Nodes appear on the element
boundaries and “fasten” the elements together.
Step 1 – Discretize the continuum
Step 2 – Select interpolation functions
Step 3 – Find the element properties
Step 4 – Assemble the element equations
Step 5 – solve the global equation systems
Step 6 – compute additional results (stress, strain, etc.)
Derivation of FEM weak forms
FEM weak forms are done when breaking down a big component for a project into
smaller parts, or preferred as the term derivation from strong to weak forms. For a certain
statement to be derived, first, the Ordinary Differential Equation (ODE) is multiplied to a
virtual function that satisfies the initial condition. Second, there will be a shifting of
derivatives. Third, integration will be done. Fourth would be the application of a boundary
condition before arriving at the weak statement of a form.
Figure 1. Graphical representation of weak forms
Galerkin’s method
This method is crucial for the Finite Element Method as this converts a differential
equation to a discrete problem. Galerkin’s method involves doing double integrations,
apply boundary conditions for the loading, and work on the equation until the deformation
or maximum deflections are found.
Discretization by triangular elements
Discretization is the process of transferring continuous functions, models,
variables and equations into discrete counterparts. The process of dividing the body into
an equivalent number of finite elements associated with nodes is called as discretization
of an element in FEM.
Post processing
After a finite element model had been prepared and checked, as well as the
boundary conditions, this phase now comes next, where investigation of the results of the
analysis will be conducted. Errors and warnings will be checked first, then reaction loads
at restrained nodes should be summed up. This basically determines whether the solution
does not contain any mathematical errors.
Boundary conditions based from known physicist has an equivalent function that it can
support. For Dirichlet boundary conditions with FEM, it is concerned with displacement.
For Neumann, it focuses on imposed strain, stress, and external loads. For Robin
boundary conditions, it can work with projects or part of projects with elastic bedding.
GIVEN ACTIVITY:
Figure 2. Activity Reference
PLOTS OF THE TWO MESHES, BOUNDARY CONDITIONS, AND LOADINGS
Figure 3. Plot for Abaqus mesh
Figure 4. Boundary Condition
PLOT THE DEFORMED SHAPE FOR BOTH MESHES ON TWO SEPARATE PLOTS
Figure 5. Figure for Plot Shape
Figure 6. Boundary Condition Structure
PLOT THE YDISPLACEMENT ALONG THE LINE ABC FOR BOTH MESHES ON THE
SAME FIGURE FOR COMPARISON
Figure 7. Fine Mesh Figure
Figure 8. MATLAB data
PLOT THE VONMISSES STRESS ON THE DEFORMED CONFIGURATION, ONLY
FOR THE FINE MESH
Figure 9. Fine Mesh Figure
Figure 10. Color Map of the Von Misses contour
TWO PRINCIPAL STRESSES FOR EACH ELEMENT IN THE COARSE MESH
Figure 11. Boundary Condition on the principal stresses
AVERAGE BAND WIDTH OF THE GLOBAL STIFFNESS MATRIX FOR BOTH
MESHES
Between the Abaqus mesh file and the MATLAB file, Abaqus was able to show the
band width effectively and the stiffness matrix through the truetoscale models and ratio,
but MATLAB was able to execute the process even though the codes are all separated
from each other, which helps easily determine which code shows a certain sector for the
program, like one for the nodal stress, displacements, elasticity graph, and the like.
PHYSICAL AND NUMERICAL IMPLICATIONS, INCLUDE CONCLUDING REMARKS
Visually and physically, the data and the models shown implied that the concrete
is elastic enough when supported by the soil, which are both able to sustain the pressure
load exerted by the water, and the dam can withstand it, up to a specific pressure load.
When viewing the model and making observations to generate data such as the stiffness
matrix, or deflections, Abaqus can effectively aid in providing the visuals. As for the
numerical implications, MATLAB is always a handy software to prove theoretical
assumptions, as well as mathematical models for topics such as this, especially for Finite
Element Analysis, which needs to focus on the data specifically as it is broken down from
the start into smaller subsets, only to be reassembled supported with equations that can
help in finding out the value of the missing variables.
BONUS QUESTION:
Implement an isoparametric triangular element and use gauss quadrature for triangles to
analyze the problem. Show that the results are in good agreement (should be identical)
with the exact integration triangle implemented in the project.
When the Gauss quadrature is integrated to the problem, it can prove that the
function over a triangular surface can be achieved by computing for the functions at
discrete points. It can be that multiplication of the discrete solutions via a weighing
function, and then running the discrete values.
When Gauss quadrature is applied to triangular areas, it states that the integral of
a function, f, over the area can be evaluated as the sum, over n integration points, of the
product of the function at each point, a weighting function for each point, and the
determinate would be:
𝑛
1
∫ 𝑓𝑑 𝑆 = ∑ 𝑤𝑖 𝑡𝑖 1𝐽𝑖
2
2
𝑖=1
APPENDIX A
The full copies of the codes are in the MATLAB and Abaqus file. For the abaqus
code, this represents what would be the neartoreal representation of the figures as well
as the mesh. The plots of the graph are represented above. As for the MATLAB code,
here we will see the actual process on how the code generates the equations and values
necessary for the FEM analysis to occur, such as the Boundary Conditions, Global
Assembly, and the Derivations of the Stiffness Matrix.
CODE FLOW:
1. Check include file, input files and adjust the values for the input:

Included_flags.m

Input_file_64ele, Input_file_16ele_tension

Input_file_16ele

Input_file_1ele
2. Preprocessor.m
3. Displacements.m, get_stress.m, gauss.m, getsctr.m, BmatElast2D, NmatElast2D
4. Mesh2D
5. Point_and_trac.m
6. Nodal_stress
7. Plotmesh
8. Solvedr.m
9. Stress_contours.m
10. Postprocessor.m
11. elasticity2d.m
Finite Element Analysis I (ENME E4332)
Fall 2020
Prof. H. Waisman
Columbia University
DEPT. OF CIVIL ENGINEERING & ENGINEERING MECHANICS
Final Project (due 12/21/2020)
Stress Analysis of a Gravity Dam on Soil Foundation
Consider the concrete built gravity dam shown in the figure. The dam is assumed to be
perfectly attached to the soil underneath. We will assume that the soil is fixed on the
bottom left end and is free to deform in the xdirection (roller boundary condition) on the
bottom right end. The dam holds back the river and hence hydrostatic pressure is
developed at the left side of the dam as shown. For simplicity the weight of the dam will be
modeled as a uniform pressure load, 𝒲, applied on its top edge.
𝒲
g
Note: all units in the figure are in meters. In addition, the following parameters are given:
Concrete
Soil
Water
Young Modulus [𝒎𝟐 ]
30 × 109
1 × 109

Poisson ratio
0.2
0.3

2200
1800
1000
𝑵
𝒌𝒈
Density [𝒎𝟑 ]
1
Finite Element Analysis I (ENME E4332)
Fall 2020
Prof. H. Waisman
Columbia University
A 2D finite element analysis is required to make sure that the dam is well designed,
assuming a plane strain constitutive law.
Assignment:
Modify the elasticity MATLAB code posted on courseworks (or write your own new code) to
analyze this problem using 3node CST triangular elements with exact integration.
The analysis should include a coarse mesh consisting of 300500 elements and a fine mesh
with 10002000 elements. Use Abaqus to generate the mesh for the model, and extract the
required element and nodal information needed to generate the same mesh in Matlab.
Notes:
1) Abaqus will generate *.inp file that will be available in your working directory, and
will include the required information needed for mesh generation.
2) When generating the mesh in Abaqus make sure you have sufficient number of nodes
placed along the line ABC where plots are required.
3) Once you are ready to submit the project, zip your MATLAB code (call the file
LastName_FirstName_UNI.zip) and upload it to Courseworks. Include a Readme file
(Readme.txt file) with the name of your team member (if any), and brief instructions
how to run the code. Your code will be tested, therefore be sure that the code can be
executed without errors before you zip it, otherwise your project will be considered
as incomplete.
4) Include with your zipped code a copy of your report in PDF or Word format. Make
sure to use legends, labels, grids and titles in all plots required for the report.
The final report (submitted as a hardcopy and also uploaded to courseworks) should
include:
1. (10%) Theoretical (physics) Background: e.g. problem statement, governing equations,
definition of plane strain and other assumptions used, anything other information you think
is important.
2
Finite Element Analysis I (ENME E4332)
Fall 2020
Prof. H. Waisman
Columbia University
2. (10%) Numerical (FEM) Background: e.g. derivation of FEM weak forms, Galerkin's
method, discretization by triangular elements, postprocessing, etc.
3. (10%) Plots of the two meshes, boundary conditions and loadings. Node numbers should
only be added to the coarse mesh figure.
4. (20%) Plot the deformed shape for both meshes on two separate plots (include the scaling
factor used for the plots – the deformation should be visible).
5. (10%) Plot the ydisplacement along the line ABC for both meshes on the same figure for
comparison.
6. (10%) Plot the VonMisses stress on the deformed configuration, only for the fine mesh.
Show the points of highest/lowest stresses. Why are these points important?
7. (10%) Compute and plot the two principal stresses for each element in the coarse mesh
(plot them as a vector field using arrows plot). Make sure to normalize the arrows in this plot,
and in addition paint (or mark) the elements with the maximum principle values.
8. (5%) The code should also output the average band width of the global stiffness matrix
(for both meshes). Add it to your report. Is the mesh numbering/connectivity ideal or can
one improve it?
9. (5%) Add a short discussion on physical and numerical implications and include some
concluding remarks.
10. (10%) Appendix A: add a brief description of the code structure, code flow and its
subroutines. The code itself should be uploaded as a zipped file to coarse works.
Bonus Question (10%):
Implement an isoparametric triangular element and use gauss quadrature for triangles to
analyze the problem. Show that the results are in good agreement (should be identical) with
the exact integration triangle implemented in the project.
Final Remark: It is strongly recommended that you prepare a nice report that may serve
you well beyond this course, for instance in job interviews.
3
...