Problems for Chapter 7
Introduction to Finite Element Analysis
in Solid Mechanics
7.2. A simple
finite element program
7.2.1. Modify the simple FEA code in FEM_conststrain.mws to
solve problems involving plane stress
deformation instead of plane strain (this should require a change to only one
line of the code). Check the modified
code by solving the problem shown in the figure. Assume that the block has unit length in both
horizontal and vertical directions, use Young’s modulus 100 and Poisson’s ratio
0.3, and take the magnitude of the distributed load to be 10 (all in arbitrary
units). Compare the predictions of the FEA analysis with the exact solution.
the simple FEA code in FEM_conststrain.mws to solve problems involving axially symmetric solids. The figure shows a representative problem to
be solved. It represents a slice through
an axially symmetric cylinder, which is prevented from stretching vertically,
and pressurized on its interior surface.
The solid is meshed using triangular elements, and the displacements are
Show that the
nonzero strain components in the element can be expressed as
Let denote the stress in the element. Find a matrix that satisfies
Write down an
expression for the strain energy density of the element.
The total strain
energy of each element must be computed.
Note that each element represents a cylindrical region of material
around the axis of symmetry. The total
strain energy in this material follows as
energy can be computed with sufficient accuracy by evaluating the integrand at
the centroid of the element, and multiplying by the area of the element, with
where denotes the radial position of the element
centroid, and is the strain energy density at the element
centroid. Use this result to deduce an
expression for the element stiffness, and modify the procedure elstif() in the
MAPLE code to compute the element stiffness.
to the potential energy from the pressure acting on element faces must also be
computed. Following the procedure
described in Chapter 7, the potential energy is
and denote the displacements at the ends of the
element face, and denote the radial position of the ends of the
element face. Calculate an expression
for P of the form
where A and B are constants that you must
determine. Modify the procedure
elresid() to implement modified element residual.
Test your routine
by calculating the stress in a pressurized cylinder, which has inner radius 1,
exterior radius 2, and is subjected to pressure p=1 on its internal bore (all in arbitrary units), and deforms
under plane strain conditions. Compare
the FEA solution for displacements and stresses with the exact solution. Run tests with different mesh densities, and
compare the results with the analytical solution.
7.2.3. Modify the simple FEA code in FEM_conststrain.mws to
solve problems which involve thermal expansion.
To this end
generic element in the mesh. Assume that
the material inside the element has a uniform thermal expansion coefficient ,
and its temperature is increased by . Let [B] and [D] denote the matrices of shape
function derivatives and material properties defined in Sections 7.2.4, and let
denote a thermal strain vector. Write down the strain energy density in the
element, in terms of these quantities and the element displacement vector .
Hence, devise a
way to calculate the total potential energy of a finite element mesh,
accounting for the effects of thermal expansion.
Modify the FEA
code to read the thermal expansion coefficient and the change in temperature
must be read from the input file, and store them as additional material
Modify the FEA
code to add the terms associated with thermal expansion to the system of
equations. It is best to do this by
writing a procedure that computes the contribution to the equation system from
one element, and then add a section to the main analysis procedure to assemble
the contributions from all elements into the global system of equations.
Test your code
using the simple test problem
7.2.4. Modify the simple FEA code in FEM_conststrain.mws to
solve plane stress problems using rectangular elements. Use the following procedure. To keep things simple, assume that the sides
of each element are parallel to the and axes, as shown in the picture. Let ,
denote the components of displacement at nodes
a, b, c, d. The displacement at an arbitrary point within
the element can be interpolated between values at the corners, as follows
that the components of nonzero
infinitesimal strain at an arbitrary point within the element may be expressed
the section of the code which reads the element connectivity, to read an extra
node for each element. To do this, you
will need to increase the size of the array named connect from
connect(1..nelem,1..3) to connect(1..nelem,1..4), and read an extra integer
node number for each element.
the procedure named elstif, which defines the element stiffness, you will need
to make the following changes. (a) You will need to modify the [B] matrix to
look like the one in Problem 11.1. Don’t forget to change the size of the [B] array from
bmat:=array(1..3,1..6) to bmat:=array(1..3,1..8). Also, note that as long as you calculate the
lengths of the element sides B and H correctly, you can use the [B] matrix
given above even if node a does not
coincide with the origin. This is
because the element stiffness only depends on the shape of the element, not on
its position. (b) To evaluate the
element stiffness, you cannot assume that is constant within the element, so instead of
multiplying by the element area, you will need to
integrate over the area of the element
Note that MAPLE will not automatically integrate each
term in a matrix. There are various ways
to fix this. One approach is to integrate each term in the matrix
then for i=1..8, j=1..8 let
Use two nested int() statements to do the
integrals. Note also that to correctly
return a matrix value for elstif, the last line of the procedure must read elstif=k,
where k is the fully assembled
Just before the
call to the elstif procedure, you will need to change the dimensions of the
element stiffness matrix from k:=array(1..6,1..6) to k:=array(1..8,1..8).
You will need to
modify the loop that assembles the global stiffness matrix to include the
fourth node in each element. To do this,
you only need to change the lines that read
for i from 1 to 3 do
for i from 1 to 4 do
the same for the j loop.
You will need to
modify the part of the routine that calculates the residual forces. The only change required is to replace the
will need to modify the procedure that calculates element strains. Now that the strains vary within the element,
you need to decide where to calculate the strains. The normal procedure would be to calculate
strains at each integration point within the element, but we used MAPLE to
evaluate the integrals when assembling the stiffness matrix, so we didn’t
define any numerical integration points.
So, in this case, just calculate the strains at the center of the
test your routine, solve the problem shown in the figure (dimensions and
material properties are in arbitrary units).
7.2.5. In this problem
you will develop and apply a finite element method to calculate the shape of a
tensioned, inextensible cable subjected to transverse loading (e.g. gravity or
wind loading). The cable is pinned at A,
and passes over a frictionless pulley at B.
A tension T is applied to the end of the cable as shown. A
(nonuniform) distributed load q(x) causes the cable to deflect by a
distance w(x) as shown. For w<<L, the potential
energy of the system may be approximated as
To develop a finite element scheme to calculate w,
divide the cable into a series of 1-D finite elements as shown. Consider a generic element of length l with nodes a, b at its ends. Assume
that the load q is uniform over the element, and assume that w
varies linearly between values ,
at the two nodes.
Write down an
expression for w at an arbitrary distance s from node a,
in terms of ,
s and l.
expression for within the element, in terms of ,
an expression for the contribution to the potential energy arising from the
element shown, and show that element contribution to the potential energy may
be expressed as
expressions for the element stiffness matrix and residual vector.
finite element mesh shown in the figure.
The loading is uniform, and each element has the same
length. The cable tension is T. Calculate the global stiffness matrix and
residual vectors for the mesh, in terms of T, L, and .
Show how the
global stiffness matrix and residual vectors must be modified to enforce the
values of w at the two intermediate nodes.