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.

 

 

 

7.2.2.       Modify 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 interpolated as

$u_i(x_1,x_2)=u_i^{(a)}{N}_a(x_1,x_2)+u_i^{(b)}{N}_b(x_1,x_2)+u_i^{(c)}{N}_c(x_1,x_2)$

where

$\begin{array}{l}{N}_a(x_1,x_2)=\frac{\left(x_2-x_2^{(b)}\right)\left(x_1^{(c)}-x_1^{(b)}\right)-\left(x_1-x_1^{(b)}\right)\left(x_2^{(c)}-x_2^{(b)}\right)}{\left(x_2^{(a)}-x_2^{(b)}\right)\left(x_1^{(c)}-x_1^{(b)}\right)-\left(x_1^{(a)}-x_1^{(b)}\right)\left(x_2^{(c)}-x_2^{(b)}\right)}\\ N_b(x_1,x_2)=\frac{\left(x_2-x_2^{(c)}\right)\left(x_1^{(a)}-x_1^{(c)}\right)-\left(x_1-x_1^{(c)}\right)\left(x_2^{(a)}-x_2^{(c)}\right)}{\left(x_2^{(b)}-x_2^{(c)}\right)\left(x_1^{(a)}-x_1^{(c)}\right)-\left(x_1^{(b)}-x_1^{(c)}\right)\left(x_2^{(a)}-x_2^{(c)}\right)}\\ N_c(x_1,x_2)=\frac{\left(x_2-x_2^{(a)}\right)\left(x_1^{(b)}-x_1^{(a)}\right)-\left(x_1-x_1^{(a)}\right)\left(x_2^{(b)}-x_2^{(a)}\right)}{\left(x_2^{(c)}-x_2^{(a)}\right)\left(x_1^{(b)}-x_1^{(a)}\right)-\left(x_1^{(c)}-x_1^{(a)}\right)\left(x_2^{(b)}-x_2^{(a)}\right)}\end{array}$

 

7.2.2.1.            Show that the nonzero strain components in the element can be expressed as

7.2.2.2.             

$\begin{array}{l}\left[\epsilon \right]=\left[\begin{array}{c}{\epsilon }_{rr}\\ {\epsilon }_{zz}\\ {\epsilon }_{\theta \theta }\\ 2{\epsilon }_{rz}\end{array}\right]\left[u\right]={\left[\begin{array}{cccccccc}{u}_r^{(a)}& u_z^{(a)}& u_r^{(b)}& u_z^{(b)}& u_r^{(c)}& u_z^{(c)}& u_r^{(d)}& u_z^{(d)}\end{array}\right]}^T\\ \left[B\right]=\left[\begin{array}{cccccc}\frac{\partial N_a}{\partial r}& 0& \frac{\partial N_b}{\partial r}& 0& \frac{\partial N_c}{\partial r}& 0\\ 0& \frac{\partial N_a}{\partial z}& 0& \frac{\partial N_b}{\partial z}& 0& \frac{\partial N_c}{\partial z}\\ \frac{N_a}{r}& 0& \frac{N_b}{r}& 0& \frac{N_c}{r}& 0\\ \frac{\partial N_a}{\partial z}& \frac{\partial N_a}{\partial r}& \frac{\partial N_b}{\partial z}& \frac{\partial N_b}{\partial r}& \frac{\partial N_c}{\partial z}& \frac{\partial N_c}{\partial r}\end{array}\right]\end{array}$

 

7.2.2.3.            Let $\sigma =\left[\begin{array}{cccc}{\sigma }_{rr}& {\sigma }_{zz}& {\sigma }_{\theta \theta }& {\sigma }_{rz}\end{array}\right]$ denote the stress in the element.  Find a matrix $[D]$ that satisfies $\left[\sigma \right]=\left[D\right]\left[\epsilon \right]$

 

7.2.2.4.            Write down an expression for the strain energy density $U^{el}$ of the element.

 

7.2.2.5.            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

$W^{el}={\displaystyle \underset{A_{el}}{\int }2\pi r{U}^{el}drdz}$

The 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 the result

$W^{el}=2\pi A_{el}\overline{r}{\overline{U}}^{el}$

where $\overline{r}$ denotes the radial position of the element centroid, and ${\overline{U}}^{el}$ 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.

 

7.2.2.6.            The contribution 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

$P=-{\displaystyle \underset{0}{\overset{L}{\int }}2\pi r{t}_i{u}_{i}ds}$

where

$u_i=u_i^{(a)}\frac{s}{L}+u_i^{(c)}\left(1-\frac{s}{L}\right)r=r^{(a)}\frac{s}{L}+r^{(c)}\left(1-\frac{s}{L}\right)$

and $u_i^{(a)},u_i^{(c)}$ denote the displacements at the ends of the element face, and $r^{(a)},r^{(c)}$ denote the radial position of the ends of the element face.  Calculate an expression for P of the form

$P^{\text{element}}=-\left[\begin{array}{cccc}{t}_{1}A& t_{2}A& t_{1}B& t_{2}B\end{array}\right]\cdot \left[\begin{array}{cccc}{u}_1^{(a)}& u_2^{(a)}& u_1^{(c)}& u_2^{(c)}\end{array}\right]$

where A and B are constants that you must determine.  Modify the procedure elresid() to implement modified element residual.

 

7.2.2.7.            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

 

7.2.3.1.            Consider a generic element in the mesh.  Assume that the material inside the element has a uniform thermal expansion coefficient $\alpha$, and its temperature is increased by $\Delta T$.  Let [B] and [D] denote the matrices of shape function derivatives and material properties defined in Sections 7.2.4, and let $\underset{\_}{q}=\alpha \Delta T{\left[1,1,0\right]}^T$ denote a thermal strain vector.  Write down the strain energy density in the element, in terms of these quantities and the element displacement vector ${\underset{\_}{u}}^{\text{element}}$.

 

7.2.3.2.            Hence, devise a way to calculate the total potential energy of a finite element mesh, accounting for the effects of thermal expansion.

 

7.2.3.3.            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 properties.

 

7.2.3.4.            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.

 

7.2.3.5.            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 $e_1$ and $e_2$ axes, as shown in the picture. Let $(u_1{}^{(a)},u_2{}^{(a)})$, $(u_1{}^{(b)},u_2{}^{(b)})$, $(u_1{}^{(c)},u_2{}^{(c)})$, $(u_1{}^{(d)},u_2{}^{(d)})$ 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

$\begin{array}{l}{u}_1=(1-\xi )(1-\eta )u_1{}^{(a)}+\xi (1-\eta )u_1{}^{(b)}+\xi \eta u_1{}^{(c)}+(1-\xi )\eta u_1{}^{(d)}\\ u_2=(1-\xi )(1-\eta )u_2{}^{(a)}+\xi (1-\eta )u_2{}^{(b)}+\xi \eta u_2{}^{(c)}+(1-\xi )\eta u_2{}^{(d)}\end{array}$

where

$\xi =x_1/B,\eta =x_2/H$

 

7.2.4.1.            Show that the components of nonzero infinitesimal strain at an arbitrary point within the element may be expressed as $\left[\epsilon \right]=[B]\left[u\right]$, where

$\begin{array}{l}\left[\epsilon \right]=\left[\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ {\epsilon }_{12}\end{array}\right]\left[u\right]={\left[\begin{array}{cccccccc}{u}_1^{(a)}& u_2^{(a)}& u_1^{(b)}& u_2^{(b)}& u_1^{(c)}& u_2^{(c)}& u_1^{(d)}& u_2^{(d)}\end{array}\right]}^T\\ \left[B\right]=\left[\begin{array}{cccccccc}-\frac{1}{B}\left(1-\frac{x_2}{H}\right)& 0& \frac{1}{B}\left(1-\frac{x_2}{H}\right)& 0& \frac{x_2}{BH}& 0& \frac{-x_2}{BH}& 0\\ 0& -\frac{1}{H}\left(1-\frac{x_1}{B}\right)& 0& \frac{-x_1}{BH}& 0& \frac{x_1}{BH}& 0& \frac{1}{H}\left(1-\frac{x_1}{B}\right)\\ -\frac{1}{2H}\left(1-\frac{x_1}{B}\right)& -\frac{1}{2B}\left(1-\frac{x_2}{H}\right)& -\frac{x_1}{2BH}& \frac{1}{2B}\left(1-\frac{x_2}{H}\right)& \frac{x_1}{2BH}& \frac{x_2}{2BH}& \frac{1}{2H}\left(1-\frac{x_1}{B}\right)& \frac{-x_2}{2BH}\end{array}\right]\end{array}$

 

7.2.4.2.            Modify 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.

 

7.2.4.3.            In 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 ${[B]}^T[D][B]$ is constant within the element, so instead of multiplying  ${[B]}^T[D][B]$ by the element area, you will need to integrate over the area of the element

$\left[K^{\text{elem}}\right]={\displaystyle \underset{0}{\overset{H}{\int }}{\displaystyle \underset{0}{\overset{B}{\int }}{\left[B\right]}^T\left[D\right]\left[B\right]}d{x}_{1}d{x}_2}$

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 separately.  Let

$[A]={\left[B\right]}^T[D]\left[B\right]$

then for i=1..8, j=1..8 let

$k_{ij}^{\text{elem}}={\displaystyle \underset{0}{\overset{H}{\int }}{\displaystyle \underset{0}{\overset{B}{\int }}{a}_{ij}}d{x}_{1}d{x}_2}$

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 stiffness matrix.

 

7.2.4.4.            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).

 

7.2.4.5.            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

to

              for i from 1 to 4  do

and the same for the j loop.

 

7.2.4.6.            You will need to modify the part of the routine that calculates the residual forces.  The only change required is to replace the line reading

>pointer := array(1..3,[2,3,1]):

with

>pointer:= array(1..4,[2,3,4,1]):

 

7.2.4.7.            You 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 element. 

 

7.2.4.8.            To 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

$V(w)={{\displaystyle \underset{0}{\overset{L}{\int }}\frac{T}{2}\left(\frac{dw}{dx}\right)}}^{2}dx-{\displaystyle \underset{0}{\overset{L}{\int }}qwdx}$

 

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 $w_a$, $w_b$ at the two nodes.

 

7.2.5.1.            Write down an expression for w at an arbitrary distance s from node a, in terms of $w_a$, $w_b$, s and l.

7.2.5.2.            Deduce an expression for $dw/dx$ within the element, in terms of $w_a$, $w_b$ and l

7.2.5.3.            Hence, calculate 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

$V^{\text{elem}}=\frac{1}{2}\left[w_a,w_b\right]\left[\begin{array}{cc}T/l& -T/l\\ -T/l& T/l\end{array}\right]\left[\begin{array}{c}{w}_a\\ w_b\end{array}\right]-\left[w_a,w_b\right]\left[\begin{array}{c}ql/2\\ ql/2\end{array}\right]$

 

7.2.5.4.            Write down expressions for the element stiffness matrix and residual vector.

 

7.2.5.5.            Consider the finite element mesh shown in the figure.  The loading $q_0$ 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 $q_0$.

 

7.2.5.6.            Show how the global stiffness matrix and residual vectors must be modified to enforce the constraints $w_1=w_4=0$

 

7.2.5.7.            Hence, calculate values of w at the two intermediate nodes.