Finite Element Analysis: Examples and Problems
Comparison of Different Elements Behaviour Under Bending
To compare the different elements described earlier, the simply supported beam with the distributed load shown in Figure 1 was modelled in the finite element analysis software ABAQUS with various different element types. The beam is supporting a distributed load and has a Young’s modulus and a length . The beam was modelled as a 2D plane shell and meshed using 2D plane stress solid elements. Plane stress assumes that the thickness of the beam is small, allowing the material to freely deform in the third direction, thereby resulting in a zero stress components in the third direction . The third direction is the 150mm dimension in this case. The thickness of the plane stress element was set to 150mm, while the value of the pressure load applied was set to . The imposed boundary conditions are at one end and a roller support at the other end. Since this is a plane problem, specifying is redundant. The beam was modelled with various types of elements and mesh sizes. The maximum normal stress components at the top and bottom fibers of the beam at midspan and the maximum vertical displacement were determined in response to the applied distributed load.
The following table compares the results for different elements with different mesh sizes measured by the number of elements (layers) in the direction of the second basis vector. Model 11 is a very fine mesh version of model 9 to show the effects of mesh refinement.
The results according to the Euler Bernoulli beam theory are as follows. Since no load is applied in the horizontal direction, the reaction in direction is equal to zero. The vertical reaction at each end can be calculated as follows:
The reaction forces in all models matched the one calculated above. The deformed shape (Figure 2) in all the models matches the shape that is expected based on the loading and boundary conditions.
The maximum displacement (obtained in the midspan) using the Euler Bernoulli beam theory for simply supported beam with a distributed load are given by:
where is the moment of inertia for the rectangular cross section give by:
The maximum bending moment for a simply supported beam is located at midspan and is given by the following equation:
The normal stress component at the bottom fibre corresponding to the maximum bending moment can be obtained as:
The Euler Bernoulli beam model predicts a linear stress profile with at the top fibres. Notice that these results are not necessarily very accurate since the Euler Bernoulli beam theory assumes that plane sections prependicular to the neutral axis before deformation remain plane and perpendicular to the neutral axis after deformation.
Behaviour of Linear Triangular Elements
The linear triangular elements dramatically underestimate the stress and the displacement. These are constant strain elements, which is evident when viewing the stress contour plot without averaging across the elements (Figure 3). Linear triangular elements are not appropriate for bending (particularly when only one layer is used) because the strain in bending is not constant, but varies linearly from the top edge to the bottom edge. It is evident from the displacement that these elements produce a very stiff structure when only one layer is used. The result is improved when three layers of elements are used because the strain is forced to be constant over a smaller area, as opposed to constant across the entire cross section of the structure. Using these elements with a very fine mesh (60 layers) comes closer to the beam theory solution with and .
Behaviour of 4Node Quadrilateral Elements
4node quadrilateral elements offer an improved solution over the linear triangular elements; however, they are still relatively stiff due to shear locking (parasitic shear) described when the element was presented here. These elements allow the stress to vary linearly within the element. In that case, they allow a negative stress at the top and a positive stress at the bottom (Figure 4). Mesh refinement to three layers greatly improves the solution, producing results that approach the Euler Bernoulli beam solution.
Effect of Using Reduced Integration
Using one layer of the 4node quadrilateral elements with reduced integration produces an extremely soft structure with wildly inaccurate results. The 4node quadrilateral elements have
4 integration points. Using reduced integration, the number of integration points is reduced to one. This is not enough to capture the actual behaviour of the element and results in a very soft structure. The integration point is at the center of the element, which is at the neutral axis of the beam when one layer of elements is used. The stresses and strains in the models with a roller support at the right end (engineering beam theory) have zero stress at the neutral axis. So, the results suggest that the elements have zero (or close to zero) stress everywhere and an extremely high displacement. It is clear that a coarse mesh of the 4node quadrilateral elements with reduced integration cannot be used to model a beam under bending. Refining the mesh to three layers produces much more reasonable results; however, the displacement is still overestimated, meaning that the modelled structure is still softer than the exact solution.
Behaviour of 8Node Quadriateral Elements
8node quadrilateral elements produce very good results, even with only one layer of elements. The higher number of nodes and integration points allows these elements to model the stress distribution within the beam with only one element in the cross section. Mesh refinement to three layers produces a slightly softer structure, with results very close to the EulerBernoulli beam solution. Using reduced integration with the 8node quadrilateral elements reduces the number of integration points from 9 to 4 with very little change in the results. Since there is very little change in the results, it would be advantageous to use reduced integration with these elements because the computational time is dramatically reduced. It is evident from the displacement results that the reduced integration 8node quadrilateral elements have essentially reached a converged solution with a threelayer mesh. When compared to a 60layer mesh (a huge increase in number of elements), very little change occurs in the results. It can also be noted that the change in the stress is much smaller than the change in the displacement between the two mesh sizes.
Symmetry and Concentrated Loads
In the cases where half the beam was modelled with a symmetry boundary condition imposed on the symmetry plane, the results are exactly the same as in the case with the same mesh and the full beam. This result is to be expected because the beam and the solution are symmetrical. The symmetry boundary condition that was imposed was to constrain the horizontal displacement along the entire symmetry plane ().
For the linear elastic material assumption, the equations of elasticity predict infinite values of the stress at the points where concentrated loads are applied. This indicates that the stress at such locations will never achieve convergence as the stress is unbounded. The boundary conditions used in this example impose a concentrated load at the corners of the beam, causing stress concentrations and a discontinuity in the deformation. The stresses further away from the concentrated load have converged, but since at the tip of the concentrated load, the predicted stresses from the elastic solution are infinite, then the finer the mesh used, the higher the values of the stress at this location. As shown in Figure 5, the stresses in the element where the concentrated reaction is applied are over 200 MPa, which is five times the maximum at the midspan of the beam.
Conclusions
Analyzing this simple beam problem highlights the importance of choosing appropriate elements, integration procedures, and mesh sizes. It was seen that lineartriangular elements are not appropriate in bending unless an extremely fine mesh is used. This is because of their constant strain/stress condition. 4node quadrilateral elements were seen to behave better than the triangular elements, but are still too stiff for this application when a coarse mesh is used. It was determined that the 8node quadrilateral elements produce very good results for this application, even when a coarse mesh is used. The different behaviour of these elements is a result of their different shape functions. The analysis emphasizes the importance of understanding the shape functions used with each element and understanding how the elements will behave in a given situation. While reduced integration can save on computational time, it must be applied carefully. With the 4node quadrilateral elements, reduced integration produced a very soft structure because the one integration point is not enough to capture the behaviour of a coarse mesh. It is evident that reduced integration should only be used with these elements if a finer mesh is used. The 8node quadrilateral elements with reduced integration seem to be the best option for this application. Using a threelayer mesh, the results are very accurate. Using reduced integration and taking advantage of symmetry saves computational time without compromising the results. In the present example, with the relatively small number of elements, the computational time is not a big factor; however, it would be more important in a larger, more complex model necessitating a finer mesh. In that case, it would be beneficial to use reduced integration and to take advantage of symmetry wherever possible.
Examples and Problems
Example 1
Find the stiffness matrix and the nodal loads due to a traction vector and a body forces vector in a plane stress element of a linear elastic small deformations material whose Young’s modulus = 1 unit and Poisson’s ratio = 0.3. Use isoparametric formulation with exact integration, full integration, and finally, with reduced integration. The geometry and loading are shown below. The thickness of the element is assumed to be equal to 1 unit.
Solution
The mapping functions between the spatial coordinate system and the element coordinate system are given by:
The gradient of the mapping is given by:
The matrix is given by:
The stiffness matrix is given by:
Where is the linear elastic isotropic plane stress constitutive relationship matrix. Exact integration produces the following stiffness matrix (units of force/length):
Setting , full integration produces the following stiffness matrix (units of force/length):
Reduced integration produces the following stiffness matrix (units of force/length):
Notice that the exact and full integration produce very similar results. The more distorted the element from a rectangle, the more the full integration technique would deviate from the exact
integration. The reducedintegration technique, however, produces numbers that highly deviate from the full integration technique. The structure is expected to be less stiff when the reduced integration technique is utilized.
Body Forces
The body forces vector in units of force/volume is given by:
The nodal forces due to the body forces vector are:
Using exact integration, the nodal forces (units of force) due to the body forces vector are:
Full integration produces the same results:
Using reduced integration, the solution is:
It is important to note that the total sum of forces is equal whether the exact, full, or reduced integration is used.
Traction Vectors
The traction vector (units of force/area) on the side connecting nodes 2 and 3 is given by:
The nodal forces due to the traction forces vector are:
The factor that transforms the integration from the spatial coordinate system of the coordinate system can be obtained as follows:
Using exact integration, full integration, and reduced integration produces the same result:
Notice that full integration is obtained using two integration points, while reduced integration is obtained using one integration point. It should be noted that the same results were obtained using the different integration techniques because the traction vector is constant. Different results would be obtained if the traction vector were not constant.
View Mathematica Code(*Using Isoparametric Formulation*) coordinates = {{0, 0}, {3, 1}, {35/10, 32/10}, {5/10, 3}}; t = 1; (*Plane Stress*) Ee = 1; Nu = 3/10; Cc = Ee/(1 + Nu)*{{1/(1  Nu), Nu/(1  Nu), 0}, {Nu/(1  Nu), 1/(1  Nu), 0}, {0, 0, 1/2}}; a = Polygon[{coordinates[[1]], coordinates[[2]], coordinates[[3]], coordinates[[4]], coordinates[[1]]}]; Graphics[a] Shapefun = Table[0, {i, 1, 4}]; Shapefun[[1]] = (1  xi) (1  eta)/4; Shapefun[[2]] = (1 + xi) (1  eta)/4; Shapefun[[3]] = (1 + xi) (1 + eta)/4; Shapefun[[4]] = (1  xi) (1 + eta)/4; Nn = Table[0, {i, 1, 2}, {j, 1, 8}]; Do[Nn[[1, 2 i  1]] = Nn[[2, 2 i]] = Shapefun[[i]], {i, 1, 4}]; x = FullSimplify[Sum[Shapefun[[i]]*coordinates[[i, 1]], {i, 1, 4}]]; y = FullSimplify[Sum[Shapefun[[i]]*coordinates[[i, 2]], {i, 1, 4}]]; J = Table[0, {i, 1, 2}, {j, 1, 2}]; J[[1, 1]] = D[x, xi]; J[[1, 2]] = D[x, eta]; J[[2, 1]] = D[y, xi]; J[[2, 2]] = D[y, eta]; JinvT = Transpose[Inverse[J]]; mat1 = {{1, 0, 0, 0}, {0, 0, 0, 1}, {0, 1, 1, 0}}; mat2 = {{JinvT[[1, 1]], JinvT[[1, 2]], 0, 0}, {JinvT[[2, 1]], JinvT[[2, 2]], 0, 0}, {0, 0, JinvT[[1, 1]], JinvT[[1, 2]]}, {0, 0, JinvT[[2, 1]], JinvT[[2, 2]]}}; mat3 = Table[0, {i, 1, 4}, {j, 1, 8}]; Do[mat3[[1, 2 i  1]] = mat3[[3, 2 i]] = D[Shapefun[[i]], xi]; mat3[[2, 2 i  1]] = mat3[[4, 2 i]] = D[Shapefun[[i]], eta], {i, 1, 4}]; B = Chop[FullSimplify[mat1.mat2.mat3]]; (*Exact integration*) Kbeforeintegration = FullSimplify[Transpose[B].Cc.B]; K1 = t*Integrate[Kbeforeintegration*Det[J], {xi, 1., 1.}, {eta, 1., 1.}]; K1 = Chop[K1]; K1 = Round[K1, 0.0001]; K1 // MatrixForm (*Full integration*) xiset = {1/Sqrt[3], 1/Sqrt[3]}; etaset = {1/Sqrt[3], 1/Sqrt[3]}; Kfull = t*FullSimplify[Sum[(Kbeforeintegration*Det[J] /. {xi > xiset[[i]], eta > etaset[[j]]}), {i, 1, 2}, {j, 1, 2}]]; Kfull = Round[Kfull, 0.0001]; Kfull // MatrixForm (*Reduced integration*) xiset = {0}; etaset = {0}; Kreduced = 4*t*FullSimplify[Sum[(Kbeforeintegration*Det[J] /. {xi > xiset[[i]], eta > etaset[[j]]}), {i, 1, 1}, {j, 1, 1}]]; Kreduced = Round[Kreduced, 0.0001]; Kreduced // MatrixForm (*Body Forces*) rb = {0, 1}; fbeforeintegration = t*Transpose[Nn].rb*Det[J]; f1 = Integrate[fbeforeintegration, {xi, 1, 1}, {eta, 1, 1}]; N[f1] // MatrixForm xiset = {1/Sqrt[3], 1/Sqrt[3]}; etaset = {1/Sqrt[3], 1/Sqrt[3]}; ffull = FullSimplify[Sum[(fbeforeintegration /. {xi > xiset[[i]], eta > etaset[[j]]}), {i, 1, 2}, {j, 1, 2}]]; N[ffull] // MatrixForm xiset = {0}; etaset = {0}; freduced = 4*FullSimplify[Sum[(fbeforeintegration /. {xi > xiset[[i]], eta > etaset[[j]]}), {i, 1, 1}, {j, 1, 1}]]; N[freduced] // MatrixForm (*Traction Forces*) dldeta = Sqrt[J[[1, 2]]^2 + J[[2, 2]]^2] /. xi > 1 tn = {1, 0}; fbeforeintegration = t*Transpose[Nn].tn*dldeta; f1 = Integrate[fbeforeintegration /. xi > 1, {eta, 1, 1}]; N[f1] // MatrixForm xiset = {1/Sqrt[3], 1/Sqrt[3]}; etaset = {1/Sqrt[3], 1/Sqrt[3]}; ffull = FullSimplify[Sum[(fbeforeintegration /. {xi > 1, eta > etaset[[j]]}), {j, 1, 2}]]; N[ffull] // MatrixForm xiset = {0}; etaset = {0}; freduced = 2*FullSimplify[Sum[(fbeforeintegration /. {xi > 1, eta > etaset[[j]]}), {j, 1, 1}]]; N[freduced] // MatrixForm
Example 2
The shown triangular element has a thickness of , , and . Find the stiffness matrices in the plane stress and plane strain conditions. Also, find the displacements corresponding to the shown loading and boundary conditions in the case of plane strain. Compare with the solution obtained using ABAQUS.
Solution
The isoparametric triangle element has the following shape functions:
The mapping functions between the spatial and element coordinate systems are given by:
The Jacobian of the transformation and its inverse are given by:
and . The matrix with:
Therefore, the matrix is given by:
The stiffness matrix is given by:
Full integration for both plane strain and plane stress conditions can be performed with one integration point with , and the associated combined weight is .
The following Mathematica code outputs the required stiffness matrices. The same stiffness matrices are obtained using ABAQUS. The columns of the stiffness matrix can be obtained using ABAQUS by restraining all the degrees of freedom to zero except for the degree of freedom corresponding to the desired column. The reactions obtained from ABAQUS are the columns of the stiffness matrix.
View Mathematica CodeClear[Ee, nu] nodes = {{0, 0}, {5, 0}, {2, 6}}; el1 = {1, 2, 3}; e1 = {nodes[[1]], nodes[[2]], nodes[[3]]}; Ee = 20000000000; nu = 0.2; Nn = {(1  xi  eta), xi, eta}; (*Plane Stress*) Cc2 = Ee/(1 + nu)*{{(1)/(1  nu), nu/(1  nu), 0}, {nu/(1  nu), 1/(1  nu), 0}, {0, 0, 1/2}}; (*Plane Strain*) Cc1 = Ee/(1 + nu)*{{(1  nu)/(1  2 nu), nu/(1  2 nu), 0}, {nu/(1  2 nu), (1  nu)/(1  2 nu), 0}, {0, 0, 1/2}}; xe1 = Sum[Nn[[i]]*e1[[i, 1]], {i, 1, 3}] ye1 = Sum[Nn[[i]]*e1[[i, 2]], {i, 1, 3}] J1 = FullSimplify[{{D[xe1, xi], D[xe1, eta]}, {D[ye1, xi], D[ye1, eta]}}]; Jin1T = Transpose[Inverse[J1]]; J1 // MatrixForm Jin1T // MatrixForm JD1 = Det[J1]; Mat1 = {{1, 0, 0, 0}, {0, 0, 0, 1}, {0, 1, 1, 0}}; Mat2e1 = Table[0, {i, 1, 4}, {j, 1, 4}]; Mat2e1[[1 ;; 2, 1 ;; 2]] = Jin1T; Mat2e1[[3 ;; 4, 3 ;; 4]] = Jin1T; Mat3 = Table[0, {i, 1, 4}, {j, 1, 6}]; Do[Mat3[[1, 2 i  1]] = Mat3[[3, 2 i]] = D[Nn[[i]], xi]; Mat3[[2, 2 i  1]] = Mat3[[4, 2 i]] = D[Nn[[i]], eta], {i, 1, 3}]; B1 = Mat1.Mat2e1.Mat3; B1 // MatrixForm Kb2 = Transpose[B1].Cc1.B1; Kb2 // MatrixForm; K1 = JD1*Transpose[B1].Cc1.B1; K1 // MatrixForm Ki1 = 1/2*K1 /. {xi > 1/3, eta > 1/3}; Ki1 // MatrixForm K2 = JD1*Transpose[B1].Cc2.B1; K2 // MatrixForm Ki2 = 1/2*K2 /. {xi > 1/3, eta > 1/3}; Ki2 // MatrixForm (*In the following, nodes 1 and 2 are restrained in the plane strain condition and a horizontal load of 50,000,000N and vertical of 50,000,000N are applied on node 3*) Kk = Ki1[[5 ;; 6, 5 ;; 6]]; Ff = Table[0, {i, 1, 2}]; Ff[[1]] = 50000000; Ff[[2]] = 50000000; Ff // MatrixForm u = Inverse[Kk].Ff; u // MatrixForm ud = {0, 0, 0, 0, u[[1]], u[[2]]}; ud // MatrixForm B1.ud // MatrixForm
Under the applied load, in the plane strain condition, the horizontal displacement and vertical displacement of the top node can be obtained by reducing the equations (eliminating the rows and columns corresponding to degrees of freedom , , , and ) as follows:
Therefore, , and . The corresponding strain in the element can be obtained as follows:
The same exact results for the three strains are obtained using ABAQUS (version 6.12)
Example 3
Consider the shown structure (Plane strain) with , under the shown concentrated load. Using two triangular elements, find the displacement of the top nodes.
Solution
Following the procedure in the previous example, element E1 has the following stiffness matrix with the corresponding degrees of freedom:
Element E2 has the following stiffness matrix with the corresponding degrees of freedom:
The global stiffness matrix is an matrix with the following entries and corresponding degrees of freedom:
By reducing the matrix (removing the rows and columns corresponding to , , , and , we are left with a matrix. The corresponding force vector is:
The corresponding displacements (in m.) are:
The following is the Mathematica code utilized. The variables mape1 and mape2 were used to map the local degrees of freedom of elements 1 and 2 respectively to the global degrees of freedom for the global matrix assembly.
View Mathematica Code
Clear[Ee, nu] nodes = {{0, 0}, {5, 0}, {2, 6}, {0, 6}}; el1 = {1, 2, 3}; el2 = {1, 3, 4}; e1 = {nodes[[1]], nodes[[2]], nodes[[3]]}; e2 = {nodes[[1]], nodes[[3]], nodes[[4]]}; mape1 = {2*el1[[1]]  1, 2*el1[[1]], 2*el1[[2]]  1, 2*el1[[2]], 2*el1[[3]]  1, 2*el1[[3]]}; mape2 = {2*el2[[1]]  1, 2*el2[[1]], 2*el2[[2]]  1, 2*el2[[2]], 2*el2[[3]]  1, 2*el2[[3]]}; mape1 // MatrixForm mape2 // MatrixForm Ee = 20000000000; nu = 0.2; Nn = {(1  xi  eta), xi, eta}; (*Plane Strain*) Cc1 = Ee/(1 + nu)*{{(1  nu)/(1  2 nu), nu/(1  2 nu), 0}, {nu/(1  2 nu), (1  nu)/(1  2 nu), 0}, {0, 0, 1/2}}; (*Element 1*) xe1 = Sum[Nn[[i]]*e1[[i, 1]], {i, 1, 3}] ye1 = Sum[Nn[[i]]*e1[[i, 2]], {i, 1, 3}] J1 = FullSimplify[{{D[xe1, xi], D[xe1, eta]}, {D[ye1, xi], D[ye1, eta]}}]; Jin1T = Transpose[Inverse[J1]]; JD1 = Det[J1]; Mat1 = {{1, 0, 0, 0}, {0, 0, 0, 1}, {0, 1, 1, 0}}; Mat2e1 = Table[0, {i, 1, 4}, {j, 1, 4}]; Mat2e1[[1 ;; 2, 1 ;; 2]] = Jin1T; Mat2e1[[3 ;; 4, 3 ;; 4]] = Jin1T; Mat3 = Table[0, {i, 1, 4}, {j, 1, 6}]; Do[Mat3[[1, 2 i  1]] = Mat3[[3, 2 i]] = D[Nn[[i]], xi]; Mat3[[2, 2 i  1]] = Mat3[[4, 2 i]] = D[Nn[[i]], eta], {i, 1, 3}]; B1 = Mat1.Mat2e1.Mat3; B1 // MatrixForm Kb2 = Transpose[B1].Cc1.B1; Kb2 // MatrixForm; K1 = JD1*Transpose[B1].Cc1.B1; K1 // MatrixForm; Ki1 = 1/2*K1 /. {xi > 1/3, eta > 1/3}; Ki1 // MatrixForm (*Element 2*) xe2 = Sum[Nn[[i]]*e2[[i, 1]], {i, 1, 3}] ye2 = Sum[Nn[[i]]*e2[[i, 2]], {i, 1, 3}] J2 = FullSimplify[{{D[xe2, xi], D[xe2, eta]}, {D[ye2, xi], D[ye2, eta]}}]; Jin2T = Transpose[Inverse[J2]]; JD2 = Det[J2]; Mat1 = {{1, 0, 0, 0}, {0, 0, 0, 1}, {0, 1, 1, 0}}; Mat2e2 = Table[0, {i, 1, 4}, {j, 1, 4}]; Mat2e2[[1 ;; 2, 1 ;; 2]] = Jin2T; Mat2e2[[3 ;; 4, 3 ;; 4]] = Jin2T; Mat3 = Table[0, {i, 1, 4}, {j, 1, 6}]; Do[Mat3[[1, 2 i  1]] = Mat3[[3, 2 i]] = D[Nn[[i]], xi]; Mat3[[2, 2 i  1]] = Mat3[[4, 2 i]] = D[Nn[[i]], eta], {i, 1, 3}]; B2 = Mat1.Mat2e2.Mat3; Kb2 = Transpose[B2].Cc1.B2; Kb2 // MatrixForm; K2 = JD2*Transpose[B2].Cc1.B2; K2 // MatrixForm; Ki2 = 1/2*K2 /. {xi > 1/3, eta > 1/3}; Ki2 // MatrixForm Ktotal = Table[0, {i, 1, 8}, {j, 1, 8}]; Do[Ktotal[[mape1[[i]], mape1[[j]]]] = Ktotal[[mape1[[i]], mape1[[j]]]] + Ki1[[i, j]], {i, 1, 6}, {j, 1, 6}]; Do[Ktotal[[mape2[[i]], mape2[[j]]]] = Ktotal[[mape2[[i]], mape2[[j]]]] + Ki2[[i, j]], {i, 1, 6}, {j, 1, 6}]; Ktotal // MatrixForm Kk = Ktotal[[5 ;; 8, 5 ;; 8]]; ff = {0, 0, 0, 50000000}; Inverse[Kk].ff // MatrixForm
ABAQUS results are exactly the same for this element. Below is the output from ABAQUS (version 6.12) showing the discontinuity in the longitudinal strain component .
Problems
 A plane element has length of 2 units aligned with the axis and height of 4 units that is aligned with the axis. If a pure rigidbody rotation around the (out of plane) axis of angle is applied to the element, find the corresponding variation in the magnitude of the small “apparent” strain components. Then, use a finite element analysis software to compare your theoretical results with the software results by assigning different values for . Comment on the results. Then, if the software used has a geometric nonlinearity option, inspect the effect of selecting this option on the magnitudes of the resulting strain components at the different chosen values for . (Hint: the rotation of the element can be applied by directly applying the appropriate displacements to the corner nodes.)
 Using any finite element analysis software of your choice, find the deflection at point A and the stress components at point B as a function of the number of elements used per the height of the beam. Use 4node quadrilateral full integration elements. Compare with the EulerBernoulli beam solution. Solve twice, once with Poisson’s ratio = 0 and another time with Poisson’s ratio = 0.3. Consider the thickness to be 1 units of length.
 Solve the previous problem using reduced integration 4node elements and reduced integration 8node elements. Compare with the results obtained in the previous problem. Comment on the convergence rate. Comment also on whether all the elements are converging to the same solution or not.
 Use a suitable quadrature to evaluate the following integrals and compare with the exact solution. Comment on the results in reference to the finite element analysis method integration scheme.
 Find the Jacobian matrix for a 4node quadrilateral isoparametric element whose coordinates are: (1,1), (3,2), (4,4), and (2,5). Assuming plane strain, unit thickness, and , find the stiffness matrix using (a) the builtin integration function in Mathematica, (b) full Gauss integration scheme, and (c) reduced Gauss integration scheme
 Following Example 1 above, use different element shapes to comment on the following.
 The effect of increasing the distortion of the element on and before integration.
 The effect of increasing the distortion of the element on the difference between and for the exact integration, full integration, and reduced integration schemes.
 If the isoparametric element is rectangular in shape but is rotated in space, what is the effect of the angle of rotation on and in the case of exact integration, full integration, and reduced integration schemes.
 The following are two main requirements for the shape functions of a 4node quadrilateral element that has a general nonrectangular shape:
 The sum of all the shape functions has to be equal to unity to ensure that rigid body motion is feasible.
 The displacement of the element side is fully determined only by the displacement of the nodes to which this side is connected in a manner that ensures element compatibility.
Does the isoparametric element formulation of a general element shape ensure the above requirements? Show that one or both of those requirements are not met if in Example 1 above either of the following two methods was used to find the shape functions:
Method 1: By assuming the following bilinear displacement across the element:Method 2: By assuming that the shape function of each node is equal to 1 at the specified node and zero on the two sides not connected to that node

The shown structure has , , , , and is subjected to the shown pressure. Assume that the structure has a unit thickness of .
 Solve the closed form solution of the differential equation of equilibrium assuming the only unknowns are the vertical displacements and the corresponding normal stresses ignoring the effect of Poisson’s ratio.
 Compare the solution using a finite element analysis software using and .
 The shown two dimensional plane strain linear elastic three node triangular element has two side lengths equal to 2m. , and . One of the sides is loaded with a trapezoidal load as shown in the figure.
 Find the matrix and the stiffness matrix.
 Find the nodal forces vector.
 If the nodal displacements of nodes 1, 2, and 3 of the element are given by (23,10), (0,30), and (20,0), respectively (units of mm), find the threedimensional strain and stress components at .
 Displacement finite element analyses are used to obtain the response of a structure to the applied external loading. For such problems, the term “linear” is used to designate “linear elements” and “linear response”. Briefly, explain the difference between “linear elements” and “linear response”. Give examples to justify your answer.
 Calculate the stiffness matrix associated with the spurious mode of the reduced integration 8 node quadratic solid element.
Please how can i get your Email
Hi there,
Thank you for the informative slides. These are very helpful.
Can you please send me solutions to the problems you have posted in these lecture notes?
Thank you.
Sorry, I don’t have typed solutions for these problems