Problems for Chapter 8
Theory and Implementation of the Finite Element Method
8.1. Generalized FEM for Static Linear Elasticity
8.1.1. Consider a one-dimensional isoparametric quadratic element, illustrated in the figure, and described in more detail in Section 8.1.5.
126.96.36.199. Suppose that the nodes have coordinates , , . Using a parametric plot, construct graphs showing the spatial variation of displacement in the element, assuming that the nodal displacements are given by (a) , (b) , (c) , (d)
188.8.131.52. Suppose that the nodes have coordinates , , . Plot graphs showing the spatial variation of displacement in the element for each of the four sets of nodal displacements given in 1.1.
8.1.2. Consider the finite element scheme to calculate displacements in an axially loaded 1D bar, described in Section 8.15. Calculate an exact analytical expression for the 3x3 stiffness matrix
for a quadratic 1-D element illustrated in the figure. Write a simple code to integrate the stiffness matrix numerically, using the procedure described in 8.1.5 and compare the result with the exact solution (for this test, choose material and geometric parameters that give ) Try integration schemes with 1, 2 and 3 integration points.
8.1.3. Consider the mapped, 8 noded isoparametric element illustrated in the figure. Write a simple program to plot a grid showing lines of and in the mapped element, as shown. Plot the grid with the following sets of nodal coordinates
· , ,
· , , ,
, , ,
Note that for the latter case, there is a region in the element where . This is unphysical. Consequently, if elements with curved sides are used in a mesh, they must be designed carefully to avoid this behavior. In addition, quadratic elements can perform poorly in large displacement analyses.
8.1.4. Set up an input file for the general 2D/3D linear elastic finite element code provided to test the 6 noded triangular elements. Run the test shown in the figure (dimensions and loading are in arbitrary units) and use a Young’s modulus and Poissons ratio . Compare the FEA solution to the exact solution.
8.1.5. Set up an input file for the general 2D/3D linear elastic finite element code provided to test the 4 noded tetrahedral elements. Run the test shown in the figure. Take the sides of the cube to have length 2 (arbitrary units), and take . Run the following boundary conditions:
184.108.40.206. at node 1, at node 2, and at nodes 3 and 4. The faces at subjected to uniform traction (arbitrary units)
220.127.116.11. at node 1, at node 5, and at nodes 3 and 6. The face at subjected to uniform traction (arbitrary units)
In each case compare the finite element solution with the exact solution (they should be equal)
8.1.6. Add lines to the 2D/3D linear elastic finite element code to compute the determinant and the eigenvalues of the global stiffness matrix. Calculate the determinant and eigenvalues for a mesh containing a single 4 noded quadrilateral element, with each of the boundary conditions shown in the figure. Briefly discuss the implications of the results on the nature of solutions to the finite element equations.
8.1.7. Set up an input file for the mesh shown in the figure. Use material properties and assume plane strain deformation. Run the following tests
18.104.22.168. Calculate the determinant of the global stiffness matrix
22.214.171.124. Change the code so that the element stiffness matrix is computed using only a single integration point. Calculate the determinant and eigenvalues of the global stiffness matrix.
8.1.8. Extend the 2D/3D linear elastic finite element code provided to solve problems involving anistropic elastic solids with cubic symmetry (see Section 3.2.16 for the constitutive law). This will require the following steps:
126.96.36.199. The elastic constants for the cubic crystal must be read from the input file
188.8.131.52. The orientation of the crystal must be read from the input file. The orientation of the crystal can be specified by the components of vectors parallel to the  and  crystallographic directions.
184.108.40.206. The parts of the code that compute the element stiffness matrix and stress (for post-processing) will need to be modified to use elastic constants for a cubic crystal. The calculation is complicated by the fact that the components of must be expressed in the global coordinate system, instead of a coordinate system aligned with the crystallographic directions. You will need to use the unit vectors in 6.2 to calculate the transformation matrix for the basis change, and use the basis change formulas in Section 3.2.11 to calculate .
220.127.116.11. Test your code by using it to compute the stresses and strains in a uniaxial tensile specimen made from a cubic crystal, in which the unit vectors parallel to  and  directions have components . Mesh the specimen with a single cubic 8 noded brick element, with side length 0.01m. Apply displacement boundary conditions to one face, and traction boundary conditions on another corresponding to uniaxial tension of 100MPa parallel to the axis. Run the following tests: (i) Verify that if are given values that represent an istropic material, the stresses and strains the element are independent of . (ii) Use values for representing copper (see Section 3.1.17). Define an apparent axial Young’s modulus for the specimen as . Plot a graph showing as a function of .
8.1.9. Extend the 2D/3D linear elastic finite element code to solve problems involving body forces. This will require the following steps
18.104.22.168. You will need to read a list of elements subjected to body forces, and the body force vector for each element. The list can be added to the end of the input file.
22.214.171.124. You will need to write a procedure to calculate the contribution from an individual element to the global system of finite element equations. This means evaluating the integral
over the volume of the element. The integral should be evaluated using numerical quadrature the procedure (other than the integrand) is essentially identical to computing the stiffness matrix.
126.96.36.199. You will need to add the contribution from each element to the global force vector. It is simplest to do this by modifying the procedure called globaltraction().
188.8.131.52. Test your code by using it to calculate the stress distribution in a 1D bar which is constrained as shown in the figure, and subjected to a uniform body force. Use 4 noded plane strain quadrilateral elements, and set (arbitrary units) and take the body force to have magnitude 5. Compare the displacements and stresses predicted by the finite element computation with the exact solution.
8.1.10. Implement the linear 3D wedge-shaped element shown in the figure the 2D/3D linear elastic finite element code. To construct the shape functions for the element, use the shape functions for a linear triangle to write down the variation with , and use a linear variation with . Assume that . You will need to add the shape functions and their derivatives to the appropriate procedures in the code. In addition, you will need to modify the procedures that compute the traction vectors associated with pressures acting on the element faces. Test your code by meshing a cube with two wedge-shaped elements as shown in the figure. Take the sides of the cube to have length 2 (arbitrary units), and take . Run the following boundary conditions:
184.108.40.206. at node 1, at node 2, and at nodes 3 and 4. The faces at subjected to uniform traction (arbitrary units)
220.127.116.11. at node 1, at node 5, and at nodes 3 and 6. The face at subjected to uniform traction (arbitrary units)
In each case compare the finite element solution with the exact solution (they should be equal).
(c) A.F. Bower, 2008