Beams, plates, and shafts are fundamental component parts of machines and structures often encountered in daily life. These components must function properly and effectively, resisting mechanical loadings. The integrity of these parts must be ensured by adequate strength design to avoid excessive loads and deformation. The mechanics of materials can evaluate stresses and deformations induced in components of interest, but can merely cover simple problems like a slender bar subjected to bending, a plate under tension, a cylindrical bar twisted by torsional moment, etc. The finite element method, however, can make analyses of stresses and deformations in machines and structures in realistic forms subjected to complicated loadings.
The object of this chapter is to show how to carry out finite-element analyses of fundamental types of component parts such as beams, plates, shafts, etc. subjected to loadings typical for each type of component part. The validity of the results of the finite-element analyses is discussed by comparing the results obtained by experiments or by the theory of elasticity.
Beam; Plate; Stepped bar; Bending; Tensile load; Torsion; Tensile stress; Shear stress; Circular hole; Elliptical hole; Fillet; Crack; Stress concentration; Stress intensity factor; Singularity of stress; Indentation; Hertz stress
Beams are important fundamental structural and/or machine elements; they are found in buildings and in bridges. Beams are also used as shafts in cars and trains (see Fig. 3.1), as wings in aircrafts, and as book shelves in bookstores. Arms and femurs of human beings and branches of trees are good examples of portions of living creatures that support their bodies like cantilever beams as illustrated in Fig. 3.2. Beams play important roles not only in inorganic but in organic structures.
The mechanics of beams is one of the most important subjects in engineering.
Perform an FEM analysis of a 2D cantilever beam shown in Fig. 3.3 and calculate the deflection of the beam at the loading point and the longitudinal stress distribution in the beam.
Before proceeding to the FEM analysis of the beam, let us review the solutions to the example problem obtained by the elementary beam theory. The maximum deflection of the beam δmax can be calculated by the following equation:
where l1 (= 80 mm) is the distance of the application point of the load from the rigid wall and l2 = l − l1.
The maximum tensile stress σmax(x) at x in the longitudinal direction appears at the upper surface of the beam in a cross section at x from the wall.
where l (= 90 mm) is the length, h (= 5 mm) is the height, b (= 10 mm) is the thickness, E (= 210 GPa) is Young's modulus, and I is the area moment of inertia of the cross section of the beam. For a beam having a rectangular cross section of a height h by a thickness b, the value of I can be calculated by the following equation:
Fig. 3.4 shows how to make structural analyses by using FEM. In this chapter, the analytical procedures will be explained following the flowchart illustrated in Fig. 3.4.
Here we shall analyse a rectangular slender beam of 5 mm (0.005 m) in height by 90 mm (0.09 m) in length by 10 mm (0.01 m) in width, as illustrated in Fig. 3.3. Fig. 3.5 shows the ANSYS Main Menu window where we can find layered command options imitating folders and files in the Microsoft Explorer folder window.
In order to prepare for creating the beam, the following operations should be made:
After carrying out the operations above, a window called Rectangle by 2 corners appears, as shown in Fig. 3.6, for the input of the geometry of a 2D rectangular beam.
The Rectangle by 2 corners window has four boxes for inputting the coordinates of the lower left corner point of the rectangular beam and the width and height of the beam to create. The following operations complete the creation procedure of the beam.
How to correct the shape of the model
When correcting the model, delete the area first, and repeat procedures (1) and (2) above. In order to delete the area, execute the following commands:
Then, the Delete Area and Below window opens and an upward arrow (↑) appears on the ANSYS Graphics window.
Next, we specify elastic constants of the beam. In the case of isotropic material, the elastic constants are Young's modulus and Poisson's ratio. This procedure can be performed any time before the solution procedure—for instance, after setting boundary conditions. If this procedure is missed, we cannot perform the solution procedure.
Then the Define Material Model Behavior window opens, as shown in Fig. 3.8.
Here we shall divide the beam area into finite elements. The procedures for finite-element discretisation are firstly to select the element type, secondly to input the element thickness and finally to divide the beam area into elements.
Then the Element Types window opens, as shown in Fig. 3.10.
The 8-node isoparametric element is a rectangular element which has four corner nodal points and four middle points, as shown in Fig. 3.13, and can realise the finite element analysis with higher accuracy than the 4-node linear rectangular element. The beam area is divided into these 8-node rectangular #82 finite elements.
Select [A] Real Constants in the ANSYS Main Menu as shown in Fig. 3.14.
The Global Element Sizes window opens, as shown in Fig. 3.19.
The Mesh Areas window opens, as shown in Fig. 3.20.
How to modify meshing
In the case of modifying meshing, delete the elements, and repeat procedures (1)–(4) above. Repeat procedures (1), (2), or (3) to modify the element type, change the plate thickness without changing the element type, or change the element size only.
In order to delete the elements, execute the following commands:
The Clear Areas window opens.
Here we shall impose constraint and loading conditions on nodes of the beam model. Display the nodes first to define the constraint and loading conditions.
The nodes are plotted in the ANSYS Graphics windows shown in Fig. 3.22.
An enlarged view of the models is often convenient when imposing constraint and loading conditions on the nodes. In order to zoom in the node display, execute the following commands:
The Pan-Zoom-Rotate window opens, as shown in Fig. 3.23.
Selection of nodes
The Apply U. ROT on Nodes window opens as shown in Fig. 3.25.
How to reselect nodes: Click [D] Reset button to clear the selection of the nodes before clicking the [C] OK button in procedure (2) above, and repeat procedures (1) and (2) above. The selection of nodes can be reset either by picking selected nodes after choosing the [E] Unpick button or clicking the right mouse button to turn the upward arrow upside down.
Imposing constraint conditions on nodes
The Apply U. ROT on Nodes window (see Fig. 3.27) opens after clicking the [C] OK button in procedure (2) in Section 3.1.
How to clear constraint conditions
The Delete Node Constrai… window opens.
Before imposing load conditions, click the Fit button in the Pan-Zoom-Rotate window (see Fig. 3.23) to get the whole view of the area and then zoom in the right end of the beam area for ease of the following operations.
Selection of the nodes
consecutively to open the Plot Numbering Controls window, as shown in Fig. 3.29.
and the Sort NODE Listing window opens (see Fig. 3.31). Select the [A] Coordinates only button and then click the [B] OK button.
to open the Apply F/M on Nodes window (see Fig. 3.33).
How to cancel the selection of the nodes of load application: Click the Reset button before clicking the OK button or click the right mouse button to change the upward arrow to the downward arrow and click the yellow frame. The yellow frame disappears and the selection of the node(s) of load application is cancelled.
Imposing load conditions on nodes
Click [A] OK in the Apply F/M on Nodes window to open another Apply F/M on Nodes window, as shown in Fig. 3.35.
How to delete load conditions: Execute the following commands:
to open the Delete F/M on Nodes window (see Fig. 3.37). Choose [A] FY or ALL in the Lab Force/moment to be deleted and click the OK button to delete the downward load applied to the #108 node.
The Solve Current Load Step and /STATUS Command windows appear as shown in Figs 3.38 and 3.39, respectively.
The Contour Nodal Solution Data window opens, as shown in Fig. 3.41.
Fig. 3.45 compares longitudinal stress distributions obtained by ANSYS with those by experiments and by elementary beam theory. The results obtained by three different methods agree well with one another. As the applied load increases, however, errors among the three groups of results become larger, especially at the clamped end. This tendency arises from the fact that the clamped condition can be hardly realised in the strict sense.
A stepped beam as shown in Fig. 3.50 can be created by adding two rectangular areas of different sizes:
Click the Reset button or click the right mouse button to change the upward arrow to the downward arrow and click the area(s) picked. The colour of the unpicked area(s) turns pink to light blue and the selection of the area(s) is cancelled.
A stepped beam with a rounded fillet, as shown in Fig. 3.51, can be created by subtracting a smaller rectangular area and a solid circle from a larger rectangular area:
to open the Solid Circular Area window. Input the coordinate of the centre ([A] WP X, [B] WP Y) and [C] Radius of the solid circle, and click the [D] OK button as shown in Fig. 3.56.
Area numbers can be displayed in the ANSYS Graphics window by the following procedure.
An elastic strip is subjected to distributed uniaxial tensile stress or negative pressure at one end and clamped at the other end.
Perform an FEM analysis of a 2D elastic strip subjected to a distributed stress in the longitudinal direction at one end and clamped at the other end, as shown in Fig. 3.61, and calculate the stress distributions along the cross sections at different distances from the loaded end in the strip.
In the procedures above, the geometry of the strip is input in millimetres. You must decide what kind of units to use in finite element analyses. When you input the geometry of a model to analyse in millimetres, for example, you must input applied loads in N (Newton) and Young's modulus in MPa, since 1 MPa is equivalent to 1 N/mm2. When you use metres and N as the units of length and load, respectively, you must input Young's modulus in Pa, since 1 Pa is equivalent to 1 N/m2. You can choose any system of unit you would like to, but your unit system must be consistent throughout the analyses.
Before proceeding to meshing, the right-end side of the strip area must be divided into two lines for imposing the triangular distribution of the applied stress or stress by executing the following commands.
Distributed load or stress can be defined by pressure on lines and the triangular distribution of applied load can be defined as the composite of two linear distributions of pressure which are antisymmetric to each other with respect to the centre line of the strip area.
Fig. 3.70 shows the variations of the longitudinal stress distribution in the cross section with the x-position of the elastic strip. At the right end of the strip, or at x = 200 mm, the distribution of the applied longitudinal stress takes the triangular shape, which is zero at the upper and lower corners and 10 MPa at the centre of the strip. The longitudinal stress distribution varies as the distance of the cross section from the right end of the strip increases, and the distribution becomes almost uniform at x = 180 mm, i.e. at one width distance from the end of the stress application. The total amount of stress in any cross section is the same, i.e. 2 kN in the strip and at any cross section at one width or larger distance from the end of the stress application, stress is uniformly distributed and the magnitude of stress becomes 5 MPa.
The above result is known as the principle of St Venant and is very useful in practice, in the design of structural components. Even if the stress distribution is very complicated at the loading points due to the complicated shape of load transfer equipment, one can assume a uniform stress distribution in the main parts of structural components or machine elements at some distance from the load transfer equipment.
An elastic plate with an elliptic hole in its centre subjected to uniform longitudinal tensile stress σ0 at one end and clamped at the other end. Perform an FEM analysis of a 2D elastic plate with an elliptic hole in its centre subjected to a uniform tensile stress σ0 in the longitudinal direction at one end and clamped at the other end, as shown in Fig. 3.71, and calculate the maximum longitudinal stress σmax in the plate. Calculate the stress concentration factor α = σmax/σ0 and observe the variation of the longitudinal stress distribution in the ligament between the foot of the hole and the edge of the plate.
Let us use a quarter model of the elastic plate with an elliptic hole as illustrated in Fig. 3.71, since the plate is symmetric about the horizontal and vertical centre lines. The quarter model can be created by a slender rectangular area from which an elliptic area is subtracted by using the Boolean operation described in Section A.1.
First, create the rectangular area by the following operation:
Then, create a circular area with a diameter of 10 mm, and reduce its diameter in the longitudinal direction to half of the original value to obtain the elliptic area. The following commands create a circular area by designating the coordinates (UX, UY) of the centre and the radius of the circular area:
In order to reduce the diameter of the circular area in the longitudinal direction to half of the original value, use the following Scale → Areas operation:
The circular area vanishes. Subtract the elliptic area from the rectangular area in a similar manner as described in Section A.1, i.e.
Due to the symmetry, the constraint conditions of the quarter plate model are UX-fixed condition on the left end and UY-fixed condition on the bottom side of the quarter plate model. Apply the constraint conditions onto the corresponding lines using the following commands:
Repeat the commands and operations (1)–(3) above for the bottom side of the model. Then, select UY in the Lab2 box and click the OK button in the Apply U. ROT on Lines window.
A uniform longitudinal stress can be defined by pressure on the right-end side of the plate model as follows:
Fig. 3.78 illustrates the boundary conditions applied to the centre notched plate model by the above operations.
In order to observe the variation of the longitudinal stress distribution in the ligament between the foot of the hole and the edge of the plate, carry out the following commands:
The Symbols window opens, as shown in Fig. 3.80.
The longitudinal stress reaches its maximum value at the foot of the hole and is decreased approaching to a constant value almost equal to the applied stress σ0 at some distance, say, about one major diameter distance, from the foot of the hole. This tendency can be explained by the principle of St Venant, as discussed in the previous section.
If an elliptical hole is placed in an infinite body and subjected to a uniform stress in a remote uniform stress of σ0, the maximum stress σmax occurs at the foot of the hole, i.e. Point B in Fig. 3.71, and is expressed by the following formula:
where a is the major radius, b is the minor radius, ρ is the local radius of curvature at the foot of the elliptic hole, and α0 is the stress concentration factor, which is defined as
The stress concentration factor α0 varies inversely proportional to the aspect ratio of elliptic hole b/a, namely the smaller the value of the aspect ratio b/a or the radius of curvature ρ becomes, the larger the value of the stress concentration factor α0 becomes.
In a finite plate, the maximum stress at the foot of the hole is increased due to the finite boundary of the plate. Fig. 3.82 shows the variation of the stress concentration factor α for elliptic holes having different aspect ratios with normalised major radius 2a/h, indicating that the stress concentration factor in a plate with finite width α is increased dramatically as the ligament between the foot of the hole and the plate edge becomes smaller.
From Fig. 3.82, the value of the stress intensity factor for the present elliptic hole is approximately 5.16, whereas Fig. 3.81 shows that the maximum value of the longitudinal stress obtained by the present FEM calculation is approximately 49.3 MPa, i.e. the value of the stress concentration factor is approximately 49.3/10 = 4.93. Hence, the relative error of the present calculation is approximately (4.93 − 5.16)/5.16 ≈ − 0.0446 = 4.46% which may be reasonably small.
An elastic plate with a crack of length 2a in its centre is subjected to uniform longitudinal tensile stress σ0 at one end and clamped at the other end, as illustrated in Fig. 3.83. Perform an FEM analysis of the 2D elastic plate having a crack of length 2a in its centre subjected to a uniform tensile stress σ0 in the longitudinal direction at one end and clamped at the other end, and calculate the value of the mode I (crack-opening mode) stress intensity factor.
Let us use a quarter model of the centre-cracked tension plate as illustrated in Fig. 3.83, since the plate is symmetric about the horizontal and vertical centre lines.
Here we use the singular element or the quarter point element, which can interpolate the stress distribution in the vicinity of the crack tip at which stress has the singularity where r is the distance from the crack tip (r/a ≪ 1). An ordinary isoparametric element which is familiar as Quad 8 node 82 has nodes at corners and also at the midpoint on each side of the element. A singular element, however, has the midpoint moved one quarter side distance from the original midpoint position to the node which is placed at the crack tip position. This is why the singular element is often called the quarter point element instead. ANSYS is equipped with a 2D triangular singular element only, not with 2D rectangular or 3D singular elements. Around the node at the crack tip, a circular area is created and is divided into a designated number of triangular singular elements. Each of these elements has its vertex placed at the crack tip position and has the quarter points on the two sides joining the vertex and the other two nodes.
In order to create the singular elements, the plate area must be created via key points set at the four corner points and at the crack tip position on the left-end side of the quarter plate area.
Then create the plate area via the five key points created by the procedures above, using the following commands:
Before meshing, the crack tip point around which the triangular singular elements will be created must be specified using the following commands:
The size of the meshes other than the singular elements and the elements adjacent to them can be controlled by the same procedures, as has been described in the previous sections of this chapter.
The meshing procedures are also the same as before.
Fig. 3.91 shows an enlarged view of the singular elements around the [A] crack tip, showing that six triangular elements are placed in a radial manner and that the size of the second row of elements is half the radius of the first row of elements, i.e. triangular singular elements.
Due to the symmetry, the constraint conditions of the quarter plate model are UX-fixed condition on the ligament region of the left end, i.e. the line between Key points #4 and #5, and UY-fixed condition on the bottom side of the quarter plate model. Apply these constraint conditions onto the corresponding lines by the following commands:
Repeat the commands and operations (1)–(3) above for the bottom side of the model. Then, select UY in the Lab2 box and click the OK button in the Apply U. ROT on Lines window.
A uniform longitudinal stress can be defined by pressure on the right-end side of the plate model, as described below:
Fig. 3.92 illustrates the boundary conditions applied to the centre-cracked tension plate model by the above operations.
The solution procedures are also the same as usual.
Fig. 3.95 is an enlarged view of the longitudinal stress distribution around the crack tip showing that very high tensile stress is induced at the crack tip, whereas compressive stress occurs around the crack surface, and the crack shape is parabolic.
Fig. 3.96 shows extrapolation of the values of the correction factor, or the nondimensional mode I (the crack-opening mode) stress intensity factor FI = KI/ to the point where r = 0, i.e. the crack tip position by the hybrid extrapolation method [1]. The plots in the right region of the figure are obtained by the following formula:
or
whereas those in the left region are obtained by the following formula:
or
where r is the distance from the crack tip, σx is the x-component of stress, σ is the applied uniform stress in the direction of the x-axis, a is the half crack length, E is Young's modulus, ux(θ = π), κ = 3 − 4ν for plane strain condition and κ = (3 − ν)/(1 + ν) for plane stress condition, and θ is the angle around the crack tip. The correction factor FI accounts for the effect of the finite boundary, or the edge effect of the plate. The value of FI can be evaluated either by the stress extrapolation method (the right region) or by the displacement extrapolation method (the left region). The hybrid extrapolation method [1] is a blend of the two types of the extrapolation methods above and can bring about better solutions. As illustrated in Fig. 3.95, the extrapolation lines obtained by the two methods agree well with each other when the gradient of the extrapolation line for the displacement method is rearranged by putting r′ = r/3. The value of FI obtained for the present CCT plate by the hybrid method is 1.0200.
Table 3.1 lists the values of FI calculated by the following four equations [2–6] for various values of nondimensional crack length λ = 2a/h:
Table 3.1
λ = 2a/h | Isida Eq. (3.7) | Feddersen Eq. (3.8) | Koiter Eq. (3.9) | Tada Eq. (3.10) |
---|---|---|---|---|
0.1 | 1.0060 | 1.0062 | 1.0048 | 1.0060 |
0.2 | 1.0246 | 1.0254 | 1.0208 | 1.0245 |
0.3 | 1.0577 | 1.0594 | 1.0510 | 1.0574 |
0.4 | 1.1094 | 1.1118 | 1.1000 | 1.1090 |
0.5 | 1.1867 | 1.1892 | 1.1757 | 1.1862 |
0.6 | 1.3033 | 1.3043 | 1.2921 | 1.3027 |
0.7 | 1.4882 | 1.4841 | 1.4779 | 1.4873 |
0.8 | 1.8160 | 1.7989 | 1.8075 | 1.8143 |
0.9 | 2.5776 | – | 2.5730 | 2.5767 |
Among these, Isida's solution is considered to be the best. The value of FI(0.2) by Isida is 1.0246. The relative error of the present result is (1.0200 − 1.0246)/1.0246 = − 0.45%, which is a reasonably good result.
An elastic cylinder with a radius of length a is pressed against a flat surface of a linearly elastic medium by a force P′. Perform an FEM analysis of an elastic cylinder of mild steel with a radius of length a pressed against an elastic flat plate of the same steel by a force P′ illustrated in Fig. 3.99 and calculate the resulting contact pressure induced in the flat plate.
Let us use a half model of the indentation problem illustrated in Fig. 3.99, since the problem is symmetric about the vertical centre line.
The elastic cylinder and the flat plate models must be combined into one model so as not to allow the rigid-body motions of the two bodies.
Select the 8-node isoparametric elements as usual.
In contact problems, two mating surfaces come in contact with each other exerting great force on each other. Contact elements must be used in contact problems to prevent penetration of one object into the other.
Here let us designate the element size by specifying number of divisions of lines picked. For this purpose, carry out the following commands:
The meshing procedures are the same as usual.
First select the lower surface of the quarter cylinder to which the contact elements are attached.
Repeat the Select → Entities … commands to select nodes on the lower surface of the cylinder.
Only nodes on the lower surface of the cylinder are plotted in the ANSYS Graphics window, as shown in Fig. 3.111.
Repeat the Select → Entities … commands again to select a smaller number of the nodes on the portion of the lower surface of the cylinder in the vicinity of the contact point.
Next, let us define contact elements to be attached to the lower surface of the cylinder by the commands as follows:
Then, perform the following commands to attach the contact elements to the lower surface of the cylinder.
Then, perform the following commands to attach the target elements to the top surface of the flat plate.
Due to the symmetry, the constraint conditions of the present model are in a UX-fixed condition on the left sides of the quarter cylinder and the half flat plate models, and a UY-fixed condition on the bottom side of the half flat plate model. Apply these constraint conditions to the corresponding lines by the following commands:
Repeat the commands and operations (1)–(3) above for the bottom side of the flat plate model. Then, select UY in the Lab2 box and click the OK button in the Apply U. ROT on Lines window.
Fig. 3.121 illustrates the boundary conditions applied to the FE model of the present contact problem by the above operations.
Fig. 3.123 is an enlarged contour of the y-component of stress around the contact point, showing that very high compressive stress is induced in the vicinity of the contact point; this is known as Hertz contact stress.
The open circular symbols in Fig. 3.124 show the plots of the y-component stress, or vertical stress, along the upper surface of the flat plate obtained by the present FE calculation. The solid line in the figure indicates the theoretical curve p(x) expressed by a parabola [10,11], i.e.
where x is the distance from the centre of contact, p(x) is the contact stress at a point x, k is the coefficient given by Eq. (3.16), and b is half the width of contact given by Eq. (3.17).
The maximum stress p0 is obtained at the centre of contact and expressed as
The length of contact is as small as 4.7 mm, so that the maximum value of the contact stress is as high as 271 MPa. The number of nodes that fall within the contact region is only 3, although the total number of nodes is approximately 15,000 and meshes are finer in the vicinity of the contact point. The three plots of the present results agree reasonably well with the theoretical curve, as shown in Fig. 3.124.