Table of Contents

The objective of this problem is to determine safe pressure of the pressure vessel. The pressure vessel will be analysed for two different materials – steel and aluminium. The problem will be solved in Matlab by developing a finite element code, and also using commercial finite element software. The problem will be analysed like an axisymmetric problem. In case of solving using Matlab, a cross-section of cylinder will be considered and solved for nodal displacements and elemental stresses. These results will be compared with finite element results. These elemental stresses will be checked against the yield strength of the material to see if the material is safe or fail. A factor of safety will also be included in the design.

A thick walled cylinder is subjected to a uniform internal pressure ‘p’. The cylinder has an internal radius of 6 in. and outside radius of 10 in. The cylinder is made up of two different materials. The goal is to find the safe pressure that the cylinder could withstand. Figure 1 shows the schematic of the problem defined.

A Finite Element Method (FEM) will be used to solve for the flexural torsional buckling load of a structure. FEM is a numerical method which is powerful for solving problems from many streams in engineering. FEM basically makes calculation easy for problems with results within acceptable range. In the case of linear systems, the finite element method requires the solution of a system of simultaneous equations rather than complicated differential equations.

The general steps for formulating the finite element solution are as follows:

Discretization is the process of modelling a body by dividing it into an equivalent body made up of smaller elements. For one-dimensional elements, each element will be connected to other elements at nodes where they share common points. After the body is divided into its elements, the element type to be used for each element must be selected. The structure decides the element type and the selection of element should be such that it behaves closely as the actual model.

Next, a displacement function is selected for each element. The most common displacement function is a polynomial function expressed in terms of the nodal unknowns. Number of dimensions of the element decides the total number of polynomial functions. There polynomial functions are used to find the displacement of the element. A one-dimensional element will have one displacement function, while two- and three-dimensional elements will have two and three displacement functions, and so on. The strain-displacement relationship and stress-strain relationship are then defined for each element. These relationships are necessary to derive the equations describing each finite element’s behaviour.

The stiffness matrix of an element can be derived by means of total potential energy. In the case of flexural-torsional buckling, an element stiffness matrix and an element geometric stiffness matrix will be derived from the energy equation.

The element stiffness matrices are then converted from the local to global coordinate system for the entire structure. The global stiffness matrices for each element are assembled to obtain the global stiffness matrix of the structure. With no boundary conditions in place the global stiffness matrix is a singular matrix. Application of boundary conditions ensures that the structure does not move as a rigid body. This is achieved by partitioning the global matrix into the free and restrained degrees of freedom. The section of the global stiffness matrix corresponding to the free degrees of freedom of the structure is used for solving the problem. For the flexural-torsional buckling problem, the partitioned global stiffness and geometric stiffness matrices are used to solve for the buckling loads.

The problem was solved as an axisymmetric problem. The axisymmetric problem deals with the analysis of structures of revolution under axisymmetric loading. A structure of revolution (SOR) is obtained by a generating cross section which rotates 360 degrees about an axis of revolution, as shown in figure below.

Such structures are said to be rotationally symmetric.

The problem was considered as internally pressurized cylinder as shown below:

The following figure shows how the element was sliced for element formation.

Figure 3: Slicing the Geometry ^{[1]}

The domain is subdivided into two elements and six nodes as shown in figure. The total force due to the distributed load is divided equally between nodes 3 and 6.

Figure 4: Discretization of Elements

The units used in the MATLAB calculations are kN and meter.

Table shows the element connectivity for this problem.

Element Number | Node i | Node j | Node m | Node n |

1 | 1 | 2 | 3 | 4 |

2 | 4 | 3 | 6 | 5 |

3 | 5 | 6 | 7 | 8 |

All the Mat lab coded are from the book ‘MATLAB Guide to Finite Elements by Peter I. Kattan’

The two element stiffness matrices k1,k2 and k3 are obtained by making calls to the

MATLAB function *BilinearQuadElementStiffness*. Each matrix has size 8 × 8.

Since the structure has six nodes, the size of the global stiffness matrix is 12 × 12.

Therefore to obtain K we first set up a zero matrix of size 12 × 12 then make two calls to the MATLAB function *BilinearQuadAssemble* since we have two elements in the structure. Each call to the function will assemble one element. The following are the MATLAB commands:

Further solving for elemental stress by applying forces:

Boundary conditions for Ansys are as shown in the figure:

The guide used for solving Ansys was from ‘UofA Ansys Tutorial’^{ [3]}

Ansys results for Steel are as follows:

As it is seen that this value is 19738.4 (psi) < 19800 (psi) hence it is safe. And the exact Mat Lab values were almost same at 19790 (psi).

From the Contour plot it is clear that the Max Stress for Aluminium is 3820.34 (psi). Also the value obtained from Mat Lab is 3850 (psi).

It is observed that both the results, Mat Lab codes and Ansys match closely to each other. The Mat Lab codes were close enough to the set safe Stress values. This is because it is pure calculation.

Ansys results had bigger deviation in results. It can be optimised by fine meshing. As the solution is always mesh sensitive.