- Email: [email protected]
A finite-volume approach to the stress analysis of pressurized axisymmetric structures
ELSEVIER
Inr. .I. Pres. I/es. & Piping 68 ( 19%) 3 I I-317 Copyright 0 1996 Elsevier Science Limited Printed in Northern Ireland. All rights reserved 0308~01h1/96/$l5.00
0308-0161(95)00070-4
A finite-volume approach to the stress analysis of pressurized axisymmetric structures M. A. Wheel Division
of Mechanics
of
Materials,
Department Glasgow
of Mechanical GI 123, U.K.
Engineering,
University
of Strathclyde,
(Received 28 July 1995;accepted 12 August 1995)
A finite-volume procedure for determining the displacementfields and elastic stressdistributions within structures that have axisymmetric geometries is presented.The production of a systemof linear algebraic equations from the equilibrium equations of the cells that are used to represent the structure is described in detail. An iterative technique is employed to solve these equationsand provide the displacementfield from which the strain and stress fields can be found. The procedure is used to model two axisymmetric test problems and the stressdistribution predicted by the formulation using cell meshesof increasing refinement is compared to analytical solutions and the variations predicted by various finite element procedures. Copyright 0 1996Elsevier ScienceLtd.
could not be used to represent complex regions and their use was restricted to modelling specialist problems such as the unsteady temperature and thermal stress fields around an arc welder traversing a rectangular plate’ and fast fracture in regular shaped test pieces.-’ Rapid crack propagation in pressurized pipelines has been modelled successfully by employing cells that are formulated in a cylindrical coordinate system,’ but the use of this formulation is restricted to cylindrical geometries because the cells are regularly shaped. Recently structural finite-volume procedures capable of analyzing more complicated geometries have been formulated.@ These formulations use shape functions analogous to those used by finite elements to describe the variation in displacements across the cell boundaries. An alternative approach, which generalized the rectangular cell-based formulation by employing coordinate transformations to produce arbitrarily shaped quadrilaterals that could represent more complex two-dimensional regions, has also been explored.’ For a test problem this formulation gave a representation of the stress distribution that was more accurate than those that were given by constant strain triangle and four-noded
1 INTRODUCTION
The finite-element method has dominated the field of numerical pressure vessel stress analysis for the past 25 years or more. However, an alternative approach based upon the finitevolume method, originally developed by the computational fluid dynamics community, is now being applied to certain specific problems in pressurized system design like the analysis of rapid crack propagation in gas distribution pipelines.’ Unlike the finite-element method, in which the relevant conservation principle, equilibrium of forces, is only satisfied in a global sense, the finite-volume procedure is conservative in that it guarantees that the cells used to discretize the region of interest are themselves in equilibrium and tractions are continuous across inter-cell boundaries. The great strength of the finite-element method is its geometric versatility; the response of complex geometries to applied loads can be modelled using a series of standard elements, whose capabilities are well known and widely documented. The finite-volume method, on the other hand, lacked this versatility. The rectangular cell-based formulations developed previously 311
312
M. A. Wheel
quadrilateral finite elements. The formulation described in this paper develops this approach for use in modelling axisymmetric geometries. The finite-volume approach has potentially important applications in pressure vessel design particularly in the analysis of problems where there may be interactions between the contained or surrounding fluid media and the containment structure. Examples of such problems include the transfer of thermal energy between fluids and structures, which occur in components such as heat exchangers, and flow-induced vibration, which may be encountered in piping systems. Another possible area of application of the finite-volume method in pressure system design is in the estimation of limit loads. Since the finite-volume method is conservative a nonlinear formulation would produce an equilibrium stress field from which the lower bound limit load could be estimated. This paper concentrates on an initial attempt to formulate an elastic axisymmetric finite-volume method capable of modelling realistic structures. The formulation is described in detail and then the precision to which it is able to predict stress distributions is compared to that of the displacement-based finite-element method and analytical solutions for two test problems. 2 AXISYMMETRIC FORMULATION
FINITE-VOLUME
2.1 Formulation of discretized equilibrium equations for internal cells Figure 1 shows an axisymmetric quadrilateral cell located within a cylindrical coordinate system. The cell, centre P, has four corner nodes numbered in a counterclockwise sense and eight neighbouring cells numbered in a similar manner as illustrated in Fig. 2. The sides of the cell have
L, I Fig. 1. An axisymmetric cell located within the cylindrical
coordinate system.
Node N3
-r Node N1
Node N4
Fig. 2. Coordinate systemdefinition and geometry specifica-
tion for face 1 of an internal cell.
lengths L, to L, and these are oriented at angles of (Y, to (Yerelative to the axial, z, direction. The radial offset from the z axis to the midpoints of the four sides are r, to r4 as illustrated. If the orthogonal stress components that act in the plane
of the cell upon
the sides are (T,,~ to Go,.+
(T,,, to uzz4and z~-,to r,;, then equilibrium axial, z, direction is described by
in the
Llr,a,,, sin a!, + L,r,2,;, cos(~, + L2r2~zz2sin cyz+ L,r,z,,, cos a2 + Llr3uzz3 sin tr, + Q3:72rz3cos (Y~ + L4r4u,,4 sin a4 + L,r,z,;, cos a4 = 0
(1)
Similarly, if the circumferential stresses that act normal to the plane of the cell upon the sides are 5881 to ~ee4 and the area of the cell is A then equilibrium in the radial, r, direction is given by L,r,a,,,cosa,
+L1r,z,,,sina, + L2r25,,2 cos a2 + L2r2zIlz sin a2 +
L3r3urr3
cos a3 + L3~~:?rz2 sin QI~
+ L4r45,,4 cos a4 + - z
(see,
+ 5’ssz
L,r,z,,, sin (u4
+ 5.983
+ &K&4> = 0
(2)
For the case of linear elastic behaviour with Young’s Modulus, E, and Poisson’s Ratio, Y, the stress components contained within eqns (1) and (2) can be rewritten as functions of the appropriate orthogonal strain components and therefore as functions of the global displacements in the radial and axial directions, u and V. For example, if the strain components present on side 1 are E,~~,&zrl and &Beland the corresponding displacement functions are (au/&-),, (&/a~),
Finite-volume
approach
and (u/r), then the radial stress on side 1 will be given by
E(l - Y) = (1 + v)(l - 2v) v + ~
l-v
-u
(3)
i Y,11
Expressions relating the other stress components in eqns (1) and (2) to displacement functions are also available. If it is assumed that the displacements vary linearly between the centres of neighbouring cells and along the sides of each cell then by using an appropriate coordinate transformation the displacement derivatives can be represented by suitable difference approximations. The derivative (&/Jr), for example, can be approximated by
(4)
where S, is the offset between the centre of cell P and the centre of its neighbour, cell 1, and p, is the angle between the radial axis and the line joining these points. The trigonometric terms in eqn (4) are produced by the transformation from the local, skewed [,q, coordinate system of the face, illustrated in Fig. 2, to the global r-z-based coordinate system. The circumferential strain on face 1, (u Ir),, can be represented by
(5) If approximations similar to eqns (4) and (5) are derived for the other displacement functions on side 1 and those on the remaining three sides of the cell then each of the stress components contained within eqns (1) and (2) can be rewritten in terms of displacements by employing relationships such as eqn (3). The format of the discretized equilibrium equations will then be Ah
+
j-
+ B&.+
c
4% +
1.3.5.7
c ,-1.3.5.7
c ALU/~ k=l.7.3.4
B,?, +
c B;,v,, X=1.1.3.4
= 0 (6)
where the subscript j identifies the four cells that
313
to stress ana1.d
have a side in common with the cell P and the subscript k distinguishes the four corner nodes of cell P. The coefficients A* and B* will depend upon the cell geometry and material properties. Nodal displacements within the discretized equilibrium equations must be eliminated. If a bilinear interpolation scheme is applied to the quadrilateral region formed by the centres of the cells that surround a given node then the nodal displacements can be rewritten as a function of the surrounding cell displacements. Consequently, each of the nodal displacements within the discretized equilibrium equations can be replaced by cell displacement functions and the equations will now take the form ,‘X ]“X Apup + c A,u, + Bpur + c B,v, = 0 (7) j-l ,=I which describes the displacement within a cell in terms of the displacements within its eight neighbours. Again, the coefficients A and B will be functions of the material properties and cell geometry. Once the discretized equilibrium equations have been obtained for all of the internal cells within the region under analysis a method of applying suitable boundary conditions is required in order to find the displacement field across the whole of the region. 2.2 Application
of the boundary
conditions
One obvious method of applying boundary conditions to an internal cell which has one or more sides adjacent to a boundary is to introduce the known quantities into the equilibrium equations. If stress boundary conditions are specified the applied stresses could be substituted into eqns 1 and 2 immediately and if displacement boundary conditions are prescribed then the known displacements could be introduced into the discretized equilibrium equations. However, the drawback of this approach is that in the case where stress boundary conditions are applied the displacements on the boundary will not be determined automatically, they will have to be calculated after the internal displacement field has been determined. An alternative approach is therefore adopted which employs specially formulated line cells to transfer the boundary conditions onto internal cells that are located next to a boundary. Unfortunately, this approach introduces addi-
M. A. Wheel
314
tional degrees of freedom into the problem but it has the advantage that displacements on boundaries that have only stress boundary conditions applied will be determined during the solution procedure. A line boundary cell, B, lying next to the internal cell, P, is shown in Fig. 3. When normal and tangential displacement boundary conditions, u,,, and ~4~~ are applied to B the displacements of B in the radial and axial directions, uB and uB, are given by uB = uN cos cxB- uT sin (Ye
(8)
and uH = uN sin (Ye+ uT cos LYE
(9) When normal and tangential stresses u,., and r, are applied to B instead then these can be related to the orthogonal stress components grrrB,CT,,~ and z~..~by o-N
=
co? aB + u,,B sin’ ff B + 2rrzB sin (YBcos (Yg
urrB
(10) and zr = -U,,B sin ffg cos ffB + (T,,B sin (YBcos (YB + r,.zB(cos2aB - sin’ (YB) (11) Each of the stress components within eqns (10) and (11) can be written in terms of the appropriate strain components and hence as functions of displacements by employing the constitutive relationships such as eqn (3). Displacement derivatives like (&L/&), can be approximated by expressions such as
displacements of the surrounding cells. The hoop strain (u/r), can be represented by
‘14 t-1 rB
c> ax
s Bl
BI
COS(~Bl
(13)
rn
-
PHI)
sin (YB1
- uNi)
L HI
LB cos(aB - PB) The nodal displacements can again be eliminated by employing a bilinear interpolation to express the nodal displacements in terms of the
14B
When used in conjunction with the constitutive relations like eqn (3) the approximations for the displacement functions given by expressions such as (12) and (13) can be used to replace the stress components in eqns (10) and (11). This will result in equations which relate the unknown displacements of B to those in its neighbours and which have a similar format to eqn (7) except that the applied stresses will appear on the right-hand side. In the case where mixed boundary conditions are imposed then one of eqns (8) or (9) must be used in conjunction with the discretized version of eqns (10) or (11) to provide a pair of eqns containing the displacements of B, UB and uR. When two line boundary cells meet and form an external corner, as illustrated in Fig. 4, the displacement derivative approximations for the two cells must be modified to account for the fact that the displacement of the corner are indeterminate: they cannot be interpolated from the surrounding cell displacements using the bilinear scheme mentioned previously. Derivatives such as (&.~/Jx),, must therefore be approximated by expressions like au = (UP - L4B1) cos LyBl -
(12)
=-
COS(aB,
-
(14)
PBI)
2.3 Solution of the discretized equilibrium equations When the discretized equilibrium equations for the line boundary cells are combined with the
Nl
Fig. 3. A line boundary cell located next to an internal cell.
Fig. 4. Two line boundary cells meeting at an external corner.
Finite-volume
approach
equations derived for the internal cells a system of linear algebraic equations containing all of the unknown displacements can be assembled. This system can be represented in matrix format by
[Albl = EFI
(15)
where the elements of the matrix A are the coefficients A and B in eqn (7) u is the displacement vector and the elements of the force vector, F, are zero for internal cells and equal to the applied boundary conditions for line boundary cells. The matrix A will have the same degree of sparsity as the stiffness matrix of an assembly of four-noded quadrilateral finite elements. The computer storage requirements of the finite-volume formulation are therefore similar to that required by the equivalent finite-element procedure. For the validation cases described in the next section the system of equations was solved using an iterative solution technique based upon the biconjugate gradient algorithm without preconditioning.’ In all cases, the algorithm converged to the required degree of accuracy in less than n iterations, where n was the number of degrees of freedom involved in the problem. Once the displacement field has been determined approximations like eqns (4) and (5) can be used to recover the strain components acting on the cell sides, from which the stress components can be found by employing the constitutive relations such as eqn (3).
3 VALIDATION OF THE FINITEVOLUME FORMULATION Two test cases were used to validate the finite-volume formulation presented here. The first problem, an internally pressurized thick walled sphere, was chosen because of the availability of an analytical solution for the stress distributions through the thickness of the vessel to which the predictions of the formulation could be compared. The second test problem, a thin walled cylindrical vessel with a hemispherical end, was selected to illustrate the geometric versatility of the finite-volume method when modelling a more practical geometry and to test how accurately the formulation is able to predict the concentrated bending stresses induced at the
.3I 5
to stress analysis
intersection hemisphere.
between
the
cylinder
and
the
3.1 The thick walled sphere A 1.Om thick steel sphere mean radius 0.25 m with an internal pressure of 100MPa was modelled by quadrilateral cells and conventional displacement-based triangular and quadrilateral elements contained within the commercial finite-element code ANSYS. In order to maintain a cell aspect ratio of approximately 1a0 and employ the same number of elements in the radial and meridional directions a 15” segment of the sphere was discretized as illustrated in Fig. 5. Initially, the coarse mesh shown in Fig. 5 containing two internal cells in the radial direction was generated. Subsequently, more refined meshes were produced by increasing the number of cells in the radial direction by 2 at each increment of refinement. Equivalent meshes of four- and eight-node quadrilateral finite elements with the same nodal locations as the finite volume meshes were also generated. Meshes of three-noded triangular elements were produced by subdividing each of the four-noded elements into two triangles. Figure 6 compares the through-thickness variations in circumferential and radial stress predicted by the coarse finite-volume mesh with the analytically determined distributions.” This figure clearly illustrates that the finite-volume method is capable of providing an accurate estimate of the stress distributions despite the lack of mesh refinement. Figure 7 shows how the error in the hoop stress at the bore, determined by the finite-volume and various finite-element procedures, varies with mesh refinement. The number of degrees of freedom in the radial direction is used as a measure of the mesh refinement rather than the
I-- 0.25 m --I b-
Fig. 5. Geometry
I.0 m -4 of the thick walled problem.
sphere
validation
M. A. Wheel
;;i 100.0 SL 50.0 w i2 5
0.0 -50.0 -100.0 -150.0
I 0.65
0.9
0.95
1
1.05
1.1
1.15
Radial Distance (m) Fig. 6. Comparisionof the hoop and radial stressvariations
in the thick walled sphere predicted by finite-volume method with analytical solution.
number of elements because this gives a true measure of the problem size and therefore the computational effort required to achieve a desired level of accuracy. It should be noted that a given finite-volume mesh will contain more degrees of freedom than the equivalent quadrilateral finite-element mesh because of the introduction of the line boundary cells at the inner and outer surfaces. Figure 7 indicates that the finite-volume method is able to provide an estimate of the maximum hoop stress that is at least as accurate as the estimates provided by the four-noded and eight-noded quadrilateral elements at all levels of mesh refinement investigated. 3.2 The cylinder
sphere intersection
A 1.0 m mean radius O-0625 m thick steel cylindrical vessel capped by a steel hemispherical end having the same radius and thickness as the cylinder was also modelled using the finite-
1.0 m
Fig.
8. Geometry of the cylinder sphere intersection problem.
volume method presented here. An internal pressure of 10 MPa was applied to the vessel. The problem is shown in Fig. 8 and since the geometry is symmetric only the part illustrated need be modelled. A sufficiently long portion of the cylinder was modelled to ensure that the effect of the bending stresses set up by the displacement mismatch at the interface was negligible at the point where the restraints were applied to the cylinder. Initially a model containing two internal cells in the radial direction was generated. Additional meshes were generated by increasing the number of cells in the through-thickness direction by two at each successive refinement. In all cases the number of cells along the axis of the cylinder and around the meridian of the hemisphere was chosen so as to keep the aspect ratio of all cells approximately equal to 1.0. Similar meshes of three-noded triangular and four-noded quadrilateral finite elements were also produced using ANSYS. Figure 9 shows how the bending component 25
r
0 2
4
6
No. of Degrees of Freedom
8
10
in Radial Direction
Fig. 7. Comparisonof errors in maximumhoop stressvalues
predicted by finite-volume method and various finiteelement procedures.
0
0.1
0.2
0.3
Distance from Intersection
0.4
0.5
(m)
Fig. 9. Comparison of the bending stressvariation in the cylinder sphere intersection problem determined bY finite-volume method with analytical solution.
Finite-volume
approach
60 P 5 :
m
:
-M-Finite
\
50
-A-3
Noded
Triangle
-E-4
Ncded
Quadrilateral
40
Ed 3E $ 30 ‘gj g zzm 20 .c 10 b I5
Volume
0
0
0 2
4 NO.
a
6
of Degrees of Freedom
10
in Radial Direction
Fig. 10. Comparison of errors in maximum bending stress values predicted by finite-volume method and various
finite-element procedures.
of the axial stress at the outer surface of the cylinder, determined by the coarse finite-volume mesh, varies on traversing along the cylinder away from the intersection. The bending component of the axial stress in the region near to the intersection was extracted by finding the mean of the difference between axial stress values at the inner and outer surfaces of the cylinder. Figure 9 compares the calculated bending stress distribution with the analytical solution”’ and clearly demonstrates that the finite-volume procedure is able to determine the form of the bending stress distribution accurately and provide a reasonable estimate of the maximum value of the bending stress even when using a coarse mesh. Finally, Fig. 10 shows how the error in the maximum bending stress, calculated using the finite-volume formulation and the finite-element procedures, varies as the meshes are refined. This figure again demonstrates that the accuracy of the finite-volume method is very similar to that of the four-node quadrilateral elements.
4 CONCLUSIONS A simple finite-volume formulation for determining stress distributions in loaded structures that have an axisymmetric geometry has been presented. As well as being mathematically straightforward the finite-volume method is also conservative in that all of the cells used to discretize the region being analyzed satisfy equilibrium. This ability to provide equilibrium
to stress analysis
317
stress fields could be exploited in applications such as the determination of lower bound limit loads. The formulation has been used to determine stress distributions in two internally pressurized structures. In each case it has been demonstrated that by mesh refinement the finite-volume method converges towards the analytical solution, and at some given level of refinement the accuracy of the method is similar to that provided by higher-order four- and eight-noded quadrilateral finite elements. Most importantly this research has produced a tool which could potentially be used with existing finite-volume computational fluid dynamics software to enable the economic analysis of problems where fluids and structures interact. This will be the goal of future work. ACKNOWLEDGEMENT The author would like to acknowledge the use of ANSYS through an educational licence. REFERENCES 1. Ivankovic, A. & Williams, J. G., The finite volume analysisof linear elastic dynamic fracture problems, in Dynamic Fracture Mechanics, ed. M. H. Aliabadi. Computational Mechanics Publications, Southampton. 1995. 2. Demirdzic, I. & Martinovic, D., Finite volume methods for therm0 elastic plastic stressanalysis,Camp. Meth. in Appl. Mech. and Eng., 109 (1993) 331-349.
3. Ivankovic, A., Demirdzic, I., Williams. J. G. & Leevers, P. S., Application of the finite volume method to the analysis of dynamic fracture problems, Int. J. of Fracture. 64 (1994)3.57-371. 4. Onate, E., Cervera, M. & Zienkiewicz, 0. C., A finite volume format for structural mechanics applications, Int. .I. Num. Meth. in Eng., 37 (1994) 181-201. 5. Demirdzic, I. & Muzaferija, S., Finite volume methods for stress analysis in complex domains, Int. J. Num. Meth. in Eng., 37 (1994) 3751-3766. 6. Bailey, C. & Cross, M., A finite volume procedure to solve elastic solid mechanics problems in three dimensionson an unstructured mesh,Int. J. Num. Meth. in Eng., 38 (1994) 3751-3766. 7. Wheel, M. A., A geometrically versatile finite volume formulation for plane elastostatic stress analysis, acceptedfor publication in J. of Strain Analysis, (1994). 8. Jennings, A. & McKeown, J., Matrix Computation. Chapter 11, SecondEd., John Wiley, Chichester, 1992. 9. Young, W., Roark’s Formulas for Stress and Strain, Chapter 12, Sixth Ed.. McGraw-Hill, New York, 1989. 10. Spence, J., Background analysis and introduction to shell theory, in Pressure Vessel Design, Concepts and Principles. Ed. J. Spence& A. S. Tooth. E. and F. N. Spon, London, 1994.
Inr. .I. Pres. I/es. & Piping 68 ( 19%) 3 I I-317 Copyright 0 1996 Elsevier Science Limited Printed in Northern Ireland. All rights reserved 0308~01h1/96/$l5.00
0308-0161(95)00070-4
A finite-volume approach to the stress analysis of pressurized axisymmetric structures M. A. Wheel Division
of Mechanics
of
Materials,
Department Glasgow
of Mechanical GI 123, U.K.
Engineering,
University
of Strathclyde,
(Received 28 July 1995;accepted 12 August 1995)
A finite-volume procedure for determining the displacementfields and elastic stressdistributions within structures that have axisymmetric geometries is presented.The production of a systemof linear algebraic equations from the equilibrium equations of the cells that are used to represent the structure is described in detail. An iterative technique is employed to solve these equationsand provide the displacementfield from which the strain and stress fields can be found. The procedure is used to model two axisymmetric test problems and the stressdistribution predicted by the formulation using cell meshesof increasing refinement is compared to analytical solutions and the variations predicted by various finite element procedures. Copyright 0 1996Elsevier ScienceLtd.
could not be used to represent complex regions and their use was restricted to modelling specialist problems such as the unsteady temperature and thermal stress fields around an arc welder traversing a rectangular plate’ and fast fracture in regular shaped test pieces.-’ Rapid crack propagation in pressurized pipelines has been modelled successfully by employing cells that are formulated in a cylindrical coordinate system,’ but the use of this formulation is restricted to cylindrical geometries because the cells are regularly shaped. Recently structural finite-volume procedures capable of analyzing more complicated geometries have been formulated.@ These formulations use shape functions analogous to those used by finite elements to describe the variation in displacements across the cell boundaries. An alternative approach, which generalized the rectangular cell-based formulation by employing coordinate transformations to produce arbitrarily shaped quadrilaterals that could represent more complex two-dimensional regions, has also been explored.’ For a test problem this formulation gave a representation of the stress distribution that was more accurate than those that were given by constant strain triangle and four-noded
1 INTRODUCTION
The finite-element method has dominated the field of numerical pressure vessel stress analysis for the past 25 years or more. However, an alternative approach based upon the finitevolume method, originally developed by the computational fluid dynamics community, is now being applied to certain specific problems in pressurized system design like the analysis of rapid crack propagation in gas distribution pipelines.’ Unlike the finite-element method, in which the relevant conservation principle, equilibrium of forces, is only satisfied in a global sense, the finite-volume procedure is conservative in that it guarantees that the cells used to discretize the region of interest are themselves in equilibrium and tractions are continuous across inter-cell boundaries. The great strength of the finite-element method is its geometric versatility; the response of complex geometries to applied loads can be modelled using a series of standard elements, whose capabilities are well known and widely documented. The finite-volume method, on the other hand, lacked this versatility. The rectangular cell-based formulations developed previously 311
312
M. A. Wheel
quadrilateral finite elements. The formulation described in this paper develops this approach for use in modelling axisymmetric geometries. The finite-volume approach has potentially important applications in pressure vessel design particularly in the analysis of problems where there may be interactions between the contained or surrounding fluid media and the containment structure. Examples of such problems include the transfer of thermal energy between fluids and structures, which occur in components such as heat exchangers, and flow-induced vibration, which may be encountered in piping systems. Another possible area of application of the finite-volume method in pressure system design is in the estimation of limit loads. Since the finite-volume method is conservative a nonlinear formulation would produce an equilibrium stress field from which the lower bound limit load could be estimated. This paper concentrates on an initial attempt to formulate an elastic axisymmetric finite-volume method capable of modelling realistic structures. The formulation is described in detail and then the precision to which it is able to predict stress distributions is compared to that of the displacement-based finite-element method and analytical solutions for two test problems. 2 AXISYMMETRIC FORMULATION
FINITE-VOLUME
2.1 Formulation of discretized equilibrium equations for internal cells Figure 1 shows an axisymmetric quadrilateral cell located within a cylindrical coordinate system. The cell, centre P, has four corner nodes numbered in a counterclockwise sense and eight neighbouring cells numbered in a similar manner as illustrated in Fig. 2. The sides of the cell have
L, I Fig. 1. An axisymmetric cell located within the cylindrical
coordinate system.
Node N3
-r Node N1
Node N4
Fig. 2. Coordinate systemdefinition and geometry specifica-
tion for face 1 of an internal cell.
lengths L, to L, and these are oriented at angles of (Y, to (Yerelative to the axial, z, direction. The radial offset from the z axis to the midpoints of the four sides are r, to r4 as illustrated. If the orthogonal stress components that act in the plane
of the cell upon
the sides are (T,,~ to Go,.+
(T,,, to uzz4and z~-,to r,;, then equilibrium axial, z, direction is described by
in the
Llr,a,,, sin a!, + L,r,2,;, cos(~, + L2r2~zz2sin cyz+ L,r,z,,, cos a2 + Llr3uzz3 sin tr, + Q3:72rz3cos (Y~ + L4r4u,,4 sin a4 + L,r,z,;, cos a4 = 0
(1)
Similarly, if the circumferential stresses that act normal to the plane of the cell upon the sides are 5881 to ~ee4 and the area of the cell is A then equilibrium in the radial, r, direction is given by L,r,a,,,cosa,
+L1r,z,,,sina, + L2r25,,2 cos a2 + L2r2zIlz sin a2 +
L3r3urr3
cos a3 + L3~~:?rz2 sin QI~
+ L4r45,,4 cos a4 + - z
(see,
+ 5’ssz
L,r,z,,, sin (u4
+ 5.983
+ &K&4> = 0
(2)
For the case of linear elastic behaviour with Young’s Modulus, E, and Poisson’s Ratio, Y, the stress components contained within eqns (1) and (2) can be rewritten as functions of the appropriate orthogonal strain components and therefore as functions of the global displacements in the radial and axial directions, u and V. For example, if the strain components present on side 1 are E,~~,&zrl and &Beland the corresponding displacement functions are (au/&-),, (&/a~),
Finite-volume
approach
and (u/r), then the radial stress on side 1 will be given by
E(l - Y) = (1 + v)(l - 2v) v + ~
l-v
-u
(3)
i Y,11
Expressions relating the other stress components in eqns (1) and (2) to displacement functions are also available. If it is assumed that the displacements vary linearly between the centres of neighbouring cells and along the sides of each cell then by using an appropriate coordinate transformation the displacement derivatives can be represented by suitable difference approximations. The derivative (&/Jr), for example, can be approximated by
(4)
where S, is the offset between the centre of cell P and the centre of its neighbour, cell 1, and p, is the angle between the radial axis and the line joining these points. The trigonometric terms in eqn (4) are produced by the transformation from the local, skewed [,q, coordinate system of the face, illustrated in Fig. 2, to the global r-z-based coordinate system. The circumferential strain on face 1, (u Ir),, can be represented by
(5) If approximations similar to eqns (4) and (5) are derived for the other displacement functions on side 1 and those on the remaining three sides of the cell then each of the stress components contained within eqns (1) and (2) can be rewritten in terms of displacements by employing relationships such as eqn (3). The format of the discretized equilibrium equations will then be Ah
+
j-
+ B&.+
c
4% +
1.3.5.7
c ,-1.3.5.7
c ALU/~ k=l.7.3.4
B,?, +
c B;,v,, X=1.1.3.4
= 0 (6)
where the subscript j identifies the four cells that
313
to stress ana1.d
have a side in common with the cell P and the subscript k distinguishes the four corner nodes of cell P. The coefficients A* and B* will depend upon the cell geometry and material properties. Nodal displacements within the discretized equilibrium equations must be eliminated. If a bilinear interpolation scheme is applied to the quadrilateral region formed by the centres of the cells that surround a given node then the nodal displacements can be rewritten as a function of the surrounding cell displacements. Consequently, each of the nodal displacements within the discretized equilibrium equations can be replaced by cell displacement functions and the equations will now take the form ,‘X ]“X Apup + c A,u, + Bpur + c B,v, = 0 (7) j-l ,=I which describes the displacement within a cell in terms of the displacements within its eight neighbours. Again, the coefficients A and B will be functions of the material properties and cell geometry. Once the discretized equilibrium equations have been obtained for all of the internal cells within the region under analysis a method of applying suitable boundary conditions is required in order to find the displacement field across the whole of the region. 2.2 Application
of the boundary
conditions
One obvious method of applying boundary conditions to an internal cell which has one or more sides adjacent to a boundary is to introduce the known quantities into the equilibrium equations. If stress boundary conditions are specified the applied stresses could be substituted into eqns 1 and 2 immediately and if displacement boundary conditions are prescribed then the known displacements could be introduced into the discretized equilibrium equations. However, the drawback of this approach is that in the case where stress boundary conditions are applied the displacements on the boundary will not be determined automatically, they will have to be calculated after the internal displacement field has been determined. An alternative approach is therefore adopted which employs specially formulated line cells to transfer the boundary conditions onto internal cells that are located next to a boundary. Unfortunately, this approach introduces addi-
M. A. Wheel
314
tional degrees of freedom into the problem but it has the advantage that displacements on boundaries that have only stress boundary conditions applied will be determined during the solution procedure. A line boundary cell, B, lying next to the internal cell, P, is shown in Fig. 3. When normal and tangential displacement boundary conditions, u,,, and ~4~~ are applied to B the displacements of B in the radial and axial directions, uB and uB, are given by uB = uN cos cxB- uT sin (Ye
(8)
and uH = uN sin (Ye+ uT cos LYE
(9) When normal and tangential stresses u,., and r, are applied to B instead then these can be related to the orthogonal stress components grrrB,CT,,~ and z~..~by o-N
=
co? aB + u,,B sin’ ff B + 2rrzB sin (YBcos (Yg
urrB
(10) and zr = -U,,B sin ffg cos ffB + (T,,B sin (YBcos (YB + r,.zB(cos2aB - sin’ (YB) (11) Each of the stress components within eqns (10) and (11) can be written in terms of the appropriate strain components and hence as functions of displacements by employing the constitutive relationships such as eqn (3). Displacement derivatives like (&L/&), can be approximated by expressions such as
displacements of the surrounding cells. The hoop strain (u/r), can be represented by
‘14 t-1 rB
c> ax
s Bl
BI
COS(~Bl
(13)
rn
-
PHI)
sin (YB1
- uNi)
L HI
LB cos(aB - PB) The nodal displacements can again be eliminated by employing a bilinear interpolation to express the nodal displacements in terms of the
14B
When used in conjunction with the constitutive relations like eqn (3) the approximations for the displacement functions given by expressions such as (12) and (13) can be used to replace the stress components in eqns (10) and (11). This will result in equations which relate the unknown displacements of B to those in its neighbours and which have a similar format to eqn (7) except that the applied stresses will appear on the right-hand side. In the case where mixed boundary conditions are imposed then one of eqns (8) or (9) must be used in conjunction with the discretized version of eqns (10) or (11) to provide a pair of eqns containing the displacements of B, UB and uR. When two line boundary cells meet and form an external corner, as illustrated in Fig. 4, the displacement derivative approximations for the two cells must be modified to account for the fact that the displacement of the corner are indeterminate: they cannot be interpolated from the surrounding cell displacements using the bilinear scheme mentioned previously. Derivatives such as (&.~/Jx),, must therefore be approximated by expressions like au = (UP - L4B1) cos LyBl -
(12)
=-
COS(aB,
-
(14)
PBI)
2.3 Solution of the discretized equilibrium equations When the discretized equilibrium equations for the line boundary cells are combined with the
Nl
Fig. 3. A line boundary cell located next to an internal cell.
Fig. 4. Two line boundary cells meeting at an external corner.
Finite-volume
approach
equations derived for the internal cells a system of linear algebraic equations containing all of the unknown displacements can be assembled. This system can be represented in matrix format by
[Albl = EFI
(15)
where the elements of the matrix A are the coefficients A and B in eqn (7) u is the displacement vector and the elements of the force vector, F, are zero for internal cells and equal to the applied boundary conditions for line boundary cells. The matrix A will have the same degree of sparsity as the stiffness matrix of an assembly of four-noded quadrilateral finite elements. The computer storage requirements of the finite-volume formulation are therefore similar to that required by the equivalent finite-element procedure. For the validation cases described in the next section the system of equations was solved using an iterative solution technique based upon the biconjugate gradient algorithm without preconditioning.’ In all cases, the algorithm converged to the required degree of accuracy in less than n iterations, where n was the number of degrees of freedom involved in the problem. Once the displacement field has been determined approximations like eqns (4) and (5) can be used to recover the strain components acting on the cell sides, from which the stress components can be found by employing the constitutive relations such as eqn (3).
3 VALIDATION OF THE FINITEVOLUME FORMULATION Two test cases were used to validate the finite-volume formulation presented here. The first problem, an internally pressurized thick walled sphere, was chosen because of the availability of an analytical solution for the stress distributions through the thickness of the vessel to which the predictions of the formulation could be compared. The second test problem, a thin walled cylindrical vessel with a hemispherical end, was selected to illustrate the geometric versatility of the finite-volume method when modelling a more practical geometry and to test how accurately the formulation is able to predict the concentrated bending stresses induced at the
.3I 5
to stress analysis
intersection hemisphere.
between
the
cylinder
and
the
3.1 The thick walled sphere A 1.Om thick steel sphere mean radius 0.25 m with an internal pressure of 100MPa was modelled by quadrilateral cells and conventional displacement-based triangular and quadrilateral elements contained within the commercial finite-element code ANSYS. In order to maintain a cell aspect ratio of approximately 1a0 and employ the same number of elements in the radial and meridional directions a 15” segment of the sphere was discretized as illustrated in Fig. 5. Initially, the coarse mesh shown in Fig. 5 containing two internal cells in the radial direction was generated. Subsequently, more refined meshes were produced by increasing the number of cells in the radial direction by 2 at each increment of refinement. Equivalent meshes of four- and eight-node quadrilateral finite elements with the same nodal locations as the finite volume meshes were also generated. Meshes of three-noded triangular elements were produced by subdividing each of the four-noded elements into two triangles. Figure 6 compares the through-thickness variations in circumferential and radial stress predicted by the coarse finite-volume mesh with the analytically determined distributions.” This figure clearly illustrates that the finite-volume method is capable of providing an accurate estimate of the stress distributions despite the lack of mesh refinement. Figure 7 shows how the error in the hoop stress at the bore, determined by the finite-volume and various finite-element procedures, varies with mesh refinement. The number of degrees of freedom in the radial direction is used as a measure of the mesh refinement rather than the
I-- 0.25 m --I b-
Fig. 5. Geometry
I.0 m -4 of the thick walled problem.
sphere
validation
M. A. Wheel
;;i 100.0 SL 50.0 w i2 5
0.0 -50.0 -100.0 -150.0
I 0.65
0.9
0.95
1
1.05
1.1
1.15
Radial Distance (m) Fig. 6. Comparisionof the hoop and radial stressvariations
in the thick walled sphere predicted by finite-volume method with analytical solution.
number of elements because this gives a true measure of the problem size and therefore the computational effort required to achieve a desired level of accuracy. It should be noted that a given finite-volume mesh will contain more degrees of freedom than the equivalent quadrilateral finite-element mesh because of the introduction of the line boundary cells at the inner and outer surfaces. Figure 7 indicates that the finite-volume method is able to provide an estimate of the maximum hoop stress that is at least as accurate as the estimates provided by the four-noded and eight-noded quadrilateral elements at all levels of mesh refinement investigated. 3.2 The cylinder
sphere intersection
A 1.0 m mean radius O-0625 m thick steel cylindrical vessel capped by a steel hemispherical end having the same radius and thickness as the cylinder was also modelled using the finite-
1.0 m
Fig.
8. Geometry of the cylinder sphere intersection problem.
volume method presented here. An internal pressure of 10 MPa was applied to the vessel. The problem is shown in Fig. 8 and since the geometry is symmetric only the part illustrated need be modelled. A sufficiently long portion of the cylinder was modelled to ensure that the effect of the bending stresses set up by the displacement mismatch at the interface was negligible at the point where the restraints were applied to the cylinder. Initially a model containing two internal cells in the radial direction was generated. Additional meshes were generated by increasing the number of cells in the through-thickness direction by two at each successive refinement. In all cases the number of cells along the axis of the cylinder and around the meridian of the hemisphere was chosen so as to keep the aspect ratio of all cells approximately equal to 1.0. Similar meshes of three-noded triangular and four-noded quadrilateral finite elements were also produced using ANSYS. Figure 9 shows how the bending component 25
r
0 2
4
6
No. of Degrees of Freedom
8
10
in Radial Direction
Fig. 7. Comparisonof errors in maximumhoop stressvalues
predicted by finite-volume method and various finiteelement procedures.
0
0.1
0.2
0.3
Distance from Intersection
0.4
0.5
(m)
Fig. 9. Comparison of the bending stressvariation in the cylinder sphere intersection problem determined bY finite-volume method with analytical solution.
Finite-volume
approach
60 P 5 :
m
:
-M-Finite
\
50
-A-3
Noded
Triangle
-E-4
Ncded
Quadrilateral
40
Ed 3E $ 30 ‘gj g zzm 20 .c 10 b I5
Volume
0
0
0 2
4 NO.
a
6
of Degrees of Freedom
10
in Radial Direction
Fig. 10. Comparison of errors in maximum bending stress values predicted by finite-volume method and various
finite-element procedures.
of the axial stress at the outer surface of the cylinder, determined by the coarse finite-volume mesh, varies on traversing along the cylinder away from the intersection. The bending component of the axial stress in the region near to the intersection was extracted by finding the mean of the difference between axial stress values at the inner and outer surfaces of the cylinder. Figure 9 compares the calculated bending stress distribution with the analytical solution”’ and clearly demonstrates that the finite-volume procedure is able to determine the form of the bending stress distribution accurately and provide a reasonable estimate of the maximum value of the bending stress even when using a coarse mesh. Finally, Fig. 10 shows how the error in the maximum bending stress, calculated using the finite-volume formulation and the finite-element procedures, varies as the meshes are refined. This figure again demonstrates that the accuracy of the finite-volume method is very similar to that of the four-node quadrilateral elements.
4 CONCLUSIONS A simple finite-volume formulation for determining stress distributions in loaded structures that have an axisymmetric geometry has been presented. As well as being mathematically straightforward the finite-volume method is also conservative in that all of the cells used to discretize the region being analyzed satisfy equilibrium. This ability to provide equilibrium
to stress analysis
317
stress fields could be exploited in applications such as the determination of lower bound limit loads. The formulation has been used to determine stress distributions in two internally pressurized structures. In each case it has been demonstrated that by mesh refinement the finite-volume method converges towards the analytical solution, and at some given level of refinement the accuracy of the method is similar to that provided by higher-order four- and eight-noded quadrilateral finite elements. Most importantly this research has produced a tool which could potentially be used with existing finite-volume computational fluid dynamics software to enable the economic analysis of problems where fluids and structures interact. This will be the goal of future work. ACKNOWLEDGEMENT The author would like to acknowledge the use of ANSYS through an educational licence. REFERENCES 1. Ivankovic, A. & Williams, J. G., The finite volume analysisof linear elastic dynamic fracture problems, in Dynamic Fracture Mechanics, ed. M. H. Aliabadi. Computational Mechanics Publications, Southampton. 1995. 2. Demirdzic, I. & Martinovic, D., Finite volume methods for therm0 elastic plastic stressanalysis,Camp. Meth. in Appl. Mech. and Eng., 109 (1993) 331-349.
3. Ivankovic, A., Demirdzic, I., Williams. J. G. & Leevers, P. S., Application of the finite volume method to the analysis of dynamic fracture problems, Int. J. of Fracture. 64 (1994)3.57-371. 4. Onate, E., Cervera, M. & Zienkiewicz, 0. C., A finite volume format for structural mechanics applications, Int. .I. Num. Meth. in Eng., 37 (1994) 181-201. 5. Demirdzic, I. & Muzaferija, S., Finite volume methods for stress analysis in complex domains, Int. J. Num. Meth. in Eng., 37 (1994) 3751-3766. 6. Bailey, C. & Cross, M., A finite volume procedure to solve elastic solid mechanics problems in three dimensionson an unstructured mesh,Int. J. Num. Meth. in Eng., 38 (1994) 3751-3766. 7. Wheel, M. A., A geometrically versatile finite volume formulation for plane elastostatic stress analysis, acceptedfor publication in J. of Strain Analysis, (1994). 8. Jennings, A. & McKeown, J., Matrix Computation. Chapter 11, SecondEd., John Wiley, Chichester, 1992. 9. Young, W., Roark’s Formulas for Stress and Strain, Chapter 12, Sixth Ed.. McGraw-Hill, New York, 1989. 10. Spence, J., Background analysis and introduction to shell theory, in Pressure Vessel Design, Concepts and Principles. Ed. J. Spence& A. S. Tooth. E. and F. N. Spon, London, 1994.