- Email: [email protected]
A FORTRAN program for computing the geometric properties of plane lamina and axi-symmetric bodies
Compurm
00417949(94)00301-7
& Smctures Vol. 54. No. 5. pp. 859-863, 1995 Copyright 0 1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0045-7949195 $9.50 + o.cil
A FORTRAN PROGRAM FOR COMPUTING THE GEOMETRIC PROPERTIES OF PLANE LAMINA AND AXI-SYMMETRIC BODIES S. Rajendran Department
of Mechanical
Engineering,
Indian
Institute
of Science,
Banglore
560 012, India
(Received 17 October 1993) Abstract-A FORTRAN volume, centre of gravity
program is presented for the computation of geometric properties, such as area, and moment of inertia, of plane lamina and axisymmetric solids. Finite element
shape functions and the Gauss-Legendre quadrature method are employed for computing the necessary integrals. The input to the program is the coordinates of a set of points on the boundary of the domain. The program
is validated
with test problems.
mechanisms, etc.
1. INTRODUffION
Computation of area, centre of gravity and moment of inertia is frequently encountered in engineering design problems. Handbooks [l-6] provide a number of useful formulae for the purpose. However, the formulae listed are usually for simpler geometries, such as circle, sector, segments, triangle, square, rectangle, etc. In the case of complex geometry, the geometry is usually split into several simpler geometric entities and the formulae are applied. This procedure is often laborious and time-consuming. In addition, there are also geometries which are not composed of such simple entities, for example, cam and gear tooth profiles and airfoils. As personal computers are commonly available in every design centre, a dedicated computer program for the purpose would be of much use. Facilities to compute some of the geometric properties are available as an integral part of commercial computer-aided design/drafting packages like AUTOCAD. However, a dedicated program has obvious advantages, such as simplicity, transparency to the user as regards the method employed, low computer memory requirement etc. Such a program is not available in the published literature. In this paper, a FORTRAN program is presented for the computation of geometric properties of plane lamina and axisymmetric bodies. The program can easily be implemented on a personal computer and used for computing properties such as area, perimeter, centre of gravity (cg.) and moments of inertia of plane lamina and volume, surface area and mass moments of inertia of axisymmetric bodies. In mechanical engineering, typical applications of the program would be to compute geometric/inertial properties of flywheels, crankshafts webs, cams, pulleys, gears, workspace of a planar/spatial
non-standard
2.
cross-sections
THE METHOD
of beams,
EMPLOYED
The boundary of the plane lamina (or that of the cross-section of the axisymmetric body) is represented by a set of line segments, each of which is defined by three nodal points, as shown in Fig. 1. The coordinates of other points in the line segments are approximated using finite element shape functions. Typically, the points K, L and M in Fig. 1 represent a line segment. The coordinates of any point P of this segment are expressed in terms of that of the nodal points K, L and M as follows:
x = i Nixi ,=1
(1)
1 Niy,,
(2)
Y =
r=,
where xi and yi (i = 1,2, 3) are the X- and y-coordinates of points K, L and M, respectively, and N,s are quadratic (parabolic) shape functions [7].
N,=+(l-5)
(3)
Nl=(l--5)(1+<)
N,=&l
(4)
+5).
(5)
The analytical expressions for computing the area, centre of gravity and moment of inertia are listed in 859
S. Rajendran the second column of Table I. The expressions are numerically evaluated using the Gauss-Legendre quadrature method [7,8,9]:
dy
where n is the number of sampling points for numerical integration, 5, is the natural coordinate of jth sampling point and W, is the corresponding weight. The numerical equivalent of the expressions for computing the area, centre of gravity and moment of inertia are listed in the third column of Table 1.
x
_
dx Fig. 1. Discretization
of the boundary.
Table Quantity
I. Analytical
expressions
Analytical
Plane lamina Area,
-
A
x-coordinate
of cg.,
y - coordinate
MI.
I,,
M.I. about
y-axis,
I,,
Numerical
expression
s
.v dx
1
of c.g., ycp
x-axis,
equivalent
1 --j-xydx A <
xcg
about
and numerical
expression
-j&b A < xy’dy
_
yx ’ dx s
Cross-M.I.,
Perimeter,
1
I,,
2, s
P
,/(dx)’
x ‘Y dl
.
+ (dy)’
s Axisymmetric body Axis of revolurion: x-axis Volume of revolution, V,, Surface area, S,, x-coordinate
Polar
of c.g., xcg
mass MI.
Mass M.I. about
about
,r-axis,
x-axis,
J,,
J,,
Axis OJ revolution: y-axis Volume of revolution, V,, Surface area, S,, y-coordinate
Polar
-
of c.g., y,,
mass M.I. about
y-axis,
J,,
Mass M.1. about x-axis, J,, The subscripts in the quantities (x),, (v),, (dx/dt ), and (dyjdt ), signify that they are computed at the Gauss-Legendre sampling point 5 = {,; abbreviatton M.I. stands for moment of inertia and c.g. for centre of gravity.
Computing geometric properties 3. COMPUTER
PROGRAM
The computer program was written in FORTRAN language in double precision arithmetic and is listed in Appendix A. The nomenclature of variables is given in the program. A pair of typical input-output data is also provided at the end of the program listing. Using Gauss-Legendre quadrature, a polynomial of degree p can be integrated exactly of sampling using a :(P + 1) number points [8,9]. The expression for the mass moments of inertia of an axisymmetric body involves a polynomial of degree 9 in 5. All other expressions involve polynomials of lower degree. Five sampling points are required for exact integration of polynomials of degree 9 and are therefore used for all the computations. The program expects the input in the following form: Density N, M l,X,,Y, 2, %>Y, 3, x37 Y3 .
. . ...
861
N - 1,x,-,,Y,-, N, XN,YN> where N is the total number of points, M is the number of points on the outer boundary, and x, and yi are the coordinates of the point i. If computations for axisymmetric bodies are required, and actual density value needs to be provided, otherwise a dummy value like ‘1.0’ may be used. The program does not expect additional information as regards numbering of line segments and their connectivity with nodes. The program assumes that nodes 1,2 and 3 are connected to first line segment, the nodes 3, 4 and 5 are connected to second line segment and so on. 4. Care is to be exercised so that the middle node of an element does not appear at the corner, as shown in Fig. 2(a) and that the line segment is not distorted much, as shown in Fig. 2(b). Under such conditions, the node numbering sequence has to be changed as shown in Fig. 2(c) and (d), respectively. Alternatively, a thumbrule to avoid such troubles is that even numbers should not occur at sharp corners or bends. 5. It should be noted that the node numbering sequence should finally end at the first point itself. Thus the coordinates of the last point, i.e.
Fig. 2. Careful positioning of the middle node of an element.
S. Rajendran Y I
rounding the opening are to be numbered clockwise direction [Fig. 3(b)].
in the
4. TEST PROBLEMS
-.
Two problems are presented here for testing the program. The first problem is a circular lamina with 1 m radius. The second problem is a square lamina of 4 m side with a central hole of 1 m radius. The nodes and numbering sequence are shown in Fig. 3(a) and (b). To compute the inertial properties of axisymmetric bodies, the density is assumed to be 7860 kg m-‘. An IBM compatible 386 PC/AT has been used for the computations reported here. The quantities computed using the program have been compared with those computed using formulae available in handbooks and the agreement has been found to be very good. The results are summarized in Tables 2 and 3.
, -x
10
8
9
5.
3x Fig. 3. Node numbering
sequence
for the test problems.
the Nth point, are the same as that of first point (Fig. 1). 6. The points of the boundary are to be numbered in the anti-clockwise direction [Fig. 3(a)]. If the geometry is multiply connected, the points sur-
CONCLUDING REiMARKS
For straight line segments, linear shape functions and thereby two nodes per element are sufficient. The program, however, uses quadratic shape functions uniformly. irrespective of whether the line segment is straight or curved. The two cases could have been handled separately in the program; this was not done as it would unnecessarily complicate the program and input data reading. In representing the circle, coordinates of 17 points (eight quadratic elements) have been found to be necessary [Fig. 3(a)] for computing geometric properties within 0.1% error (Table 2). On the contrary. computation using formulae would require only the coordinates of centre and radius. Therefore, in cases where the geometry is purely composed of simpler geometries such as straight lines, circles etc., the present
Table 2. Comparison of results for test problem no. I Quantity Plane lamina A xCR YCP I,, I,., I i’ Axisymmetric
x-axis
Axis of revolution: V,, St,
y - axis
J,,
J>.,
Formula
%Deviation
3.1392401 2.0000000 3.0000000 29.037384 13.341183
3.1415927 2.0000000 3 .ooooooo 29.059732 13.351769
-0.075 0.0 0.0 -0.077 -0.079
6.2809012 18.835441
6.2831853 18.849556
-0.036 -0.075
59.173282
59.217626
-0.075
4534483.3 118.39220 4243838. I
4538142.8 118.43525 4247236.2
-0.036 -0.081 - 0.080
39.448854 78.928132 1472648.9 3604395.4
39.478417 78.956835 1473926.7 3607241.7
-0.075 -0.036 -0.087 -0.079
body
Axis of revolution: V,, s J::
J,,
Program
Computing Table
3. Comparison
Quantity
geometric of results
properties for test problem
Program
863 no. 2
Formula
% Deviation
Plane lamina
A -%g )‘ce I \\ I I, I 2’
12.86075989 4.000000000 6.000000000 483.5364669 226.3212690 308.6582374 16.00000000
12.8584074 4.000000000 6.000000000 483.450600 226.282453 308.601776 16.00000000
0.018 0.0 0.0 0.018 0.017 0.018 0.0
484.8392256 603.1857895 155457137.7 144790960.0
484.7505372 603.1857895 155430997.3 144766387.9
0.018 0.0 0.017 0.017
323.2261504 402.1238597 52826940.96 121932882.1
323.1670242 402.1238597 52818808.62 121911853.4
0.018 0.0 0.015 0.017
Axisymmetric body
Axis
ofrevolution: x-axis
“,t S,,
J,, Jb,
A.xis v. ,
of revolution - y -uxis
will be inefficient. For such cases, a subroutine purely based on formulae would appear attractive. Such additional facilities could be easily added by the user. The present program, however, has been written for general purpose usage. It always handles the problems the same way, irrespective of whether the geometry is composed of simple geometric entities or not.
program
REFERENCES
1. E. Oberg and F. D. Jones, Machinery’s Handbook, 18th Edn. Industrial Press. New York (1968). 2. T. G. Hicks and S. D. Hicks (editors), Standard Handbook qf Engineering Calculations, 2nd Edn. McGraw-Hill, New York (1986).
3. H. A. Rothbart (editor), Mechanical Design and Systems Handbook, 2nd Edn. McGraw-Hill, New York (1986). 4. M. Kutz (editor), Mechanical Engineers Handbook. Wiley, New York (1986). 5. J. E. Shigley and C. R. Mischke (editors), Standard Handbook of Machine Design, McGraw-Hill, New York (1986). 6. W. C. Young, Roark’s Formulae for Stress and Strain, 6th Edn. McGraw-Hill, New York (1989). 7. E. Hinton and D. R. J. Owen, Finife Element Programming. Academic Press, London (1977). 8. 0. C. Zienkiewicz and K. Morgan, Finite Element and Approximation. Wiley, New York (1983). 9. 0. C. Zienkiewicz, Finite Element Method, 3rd Edn. Tata-McGraw-Hill, New Delhi (1979).
00417949(94)00301-7
& Smctures Vol. 54. No. 5. pp. 859-863, 1995 Copyright 0 1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0045-7949195 $9.50 + o.cil
A FORTRAN PROGRAM FOR COMPUTING THE GEOMETRIC PROPERTIES OF PLANE LAMINA AND AXI-SYMMETRIC BODIES S. Rajendran Department
of Mechanical
Engineering,
Indian
Institute
of Science,
Banglore
560 012, India
(Received 17 October 1993) Abstract-A FORTRAN volume, centre of gravity
program is presented for the computation of geometric properties, such as area, and moment of inertia, of plane lamina and axisymmetric solids. Finite element
shape functions and the Gauss-Legendre quadrature method are employed for computing the necessary integrals. The input to the program is the coordinates of a set of points on the boundary of the domain. The program
is validated
with test problems.
mechanisms, etc.
1. INTRODUffION
Computation of area, centre of gravity and moment of inertia is frequently encountered in engineering design problems. Handbooks [l-6] provide a number of useful formulae for the purpose. However, the formulae listed are usually for simpler geometries, such as circle, sector, segments, triangle, square, rectangle, etc. In the case of complex geometry, the geometry is usually split into several simpler geometric entities and the formulae are applied. This procedure is often laborious and time-consuming. In addition, there are also geometries which are not composed of such simple entities, for example, cam and gear tooth profiles and airfoils. As personal computers are commonly available in every design centre, a dedicated computer program for the purpose would be of much use. Facilities to compute some of the geometric properties are available as an integral part of commercial computer-aided design/drafting packages like AUTOCAD. However, a dedicated program has obvious advantages, such as simplicity, transparency to the user as regards the method employed, low computer memory requirement etc. Such a program is not available in the published literature. In this paper, a FORTRAN program is presented for the computation of geometric properties of plane lamina and axisymmetric bodies. The program can easily be implemented on a personal computer and used for computing properties such as area, perimeter, centre of gravity (cg.) and moments of inertia of plane lamina and volume, surface area and mass moments of inertia of axisymmetric bodies. In mechanical engineering, typical applications of the program would be to compute geometric/inertial properties of flywheels, crankshafts webs, cams, pulleys, gears, workspace of a planar/spatial
non-standard
2.
cross-sections
THE METHOD
of beams,
EMPLOYED
The boundary of the plane lamina (or that of the cross-section of the axisymmetric body) is represented by a set of line segments, each of which is defined by three nodal points, as shown in Fig. 1. The coordinates of other points in the line segments are approximated using finite element shape functions. Typically, the points K, L and M in Fig. 1 represent a line segment. The coordinates of any point P of this segment are expressed in terms of that of the nodal points K, L and M as follows:
x = i Nixi ,=1
(1)
1 Niy,,
(2)
Y =
r=,
where xi and yi (i = 1,2, 3) are the X- and y-coordinates of points K, L and M, respectively, and N,s are quadratic (parabolic) shape functions [7].
N,=+(l-5)
(3)
Nl=(l--5)(1+<)
N,=&l
(4)
+5).
(5)
The analytical expressions for computing the area, centre of gravity and moment of inertia are listed in 859
S. Rajendran the second column of Table I. The expressions are numerically evaluated using the Gauss-Legendre quadrature method [7,8,9]:
dy
where n is the number of sampling points for numerical integration, 5, is the natural coordinate of jth sampling point and W, is the corresponding weight. The numerical equivalent of the expressions for computing the area, centre of gravity and moment of inertia are listed in the third column of Table 1.
x
_
dx Fig. 1. Discretization
of the boundary.
Table Quantity
I. Analytical
expressions
Analytical
Plane lamina Area,
-
A
x-coordinate
of cg.,
y - coordinate
MI.
I,,
M.I. about
y-axis,
I,,
Numerical
expression
s
.v dx
1
of c.g., ycp
x-axis,
equivalent
1 --j-xydx A <
xcg
about
and numerical
expression
-j&b A < xy’dy
_
yx ’ dx s
Cross-M.I.,
Perimeter,
1
I,,
2, s
P
,/(dx)’
x ‘Y dl
.
+ (dy)’
s Axisymmetric body Axis of revolurion: x-axis Volume of revolution, V,, Surface area, S,, x-coordinate
Polar
of c.g., xcg
mass MI.
Mass M.I. about
about
,r-axis,
x-axis,
J,,
J,,
Axis OJ revolution: y-axis Volume of revolution, V,, Surface area, S,, y-coordinate
Polar
-
of c.g., y,,
mass M.I. about
y-axis,
J,,
Mass M.1. about x-axis, J,, The subscripts in the quantities (x),, (v),, (dx/dt ), and (dyjdt ), signify that they are computed at the Gauss-Legendre sampling point 5 = {,; abbreviatton M.I. stands for moment of inertia and c.g. for centre of gravity.
Computing geometric properties 3. COMPUTER
PROGRAM
The computer program was written in FORTRAN language in double precision arithmetic and is listed in Appendix A. The nomenclature of variables is given in the program. A pair of typical input-output data is also provided at the end of the program listing. Using Gauss-Legendre quadrature, a polynomial of degree p can be integrated exactly of sampling using a :(P + 1) number points [8,9]. The expression for the mass moments of inertia of an axisymmetric body involves a polynomial of degree 9 in 5. All other expressions involve polynomials of lower degree. Five sampling points are required for exact integration of polynomials of degree 9 and are therefore used for all the computations. The program expects the input in the following form: Density N, M l,X,,Y, 2, %>Y, 3, x37 Y3 .
. . ...
861
N - 1,x,-,,Y,-, N, XN,YN> where N is the total number of points, M is the number of points on the outer boundary, and x, and yi are the coordinates of the point i. If computations for axisymmetric bodies are required, and actual density value needs to be provided, otherwise a dummy value like ‘1.0’ may be used. The program does not expect additional information as regards numbering of line segments and their connectivity with nodes. The program assumes that nodes 1,2 and 3 are connected to first line segment, the nodes 3, 4 and 5 are connected to second line segment and so on. 4. Care is to be exercised so that the middle node of an element does not appear at the corner, as shown in Fig. 2(a) and that the line segment is not distorted much, as shown in Fig. 2(b). Under such conditions, the node numbering sequence has to be changed as shown in Fig. 2(c) and (d), respectively. Alternatively, a thumbrule to avoid such troubles is that even numbers should not occur at sharp corners or bends. 5. It should be noted that the node numbering sequence should finally end at the first point itself. Thus the coordinates of the last point, i.e.
Fig. 2. Careful positioning of the middle node of an element.
S. Rajendran Y I
rounding the opening are to be numbered clockwise direction [Fig. 3(b)].
in the
4. TEST PROBLEMS
-.
Two problems are presented here for testing the program. The first problem is a circular lamina with 1 m radius. The second problem is a square lamina of 4 m side with a central hole of 1 m radius. The nodes and numbering sequence are shown in Fig. 3(a) and (b). To compute the inertial properties of axisymmetric bodies, the density is assumed to be 7860 kg m-‘. An IBM compatible 386 PC/AT has been used for the computations reported here. The quantities computed using the program have been compared with those computed using formulae available in handbooks and the agreement has been found to be very good. The results are summarized in Tables 2 and 3.
, -x
10
8
9
5.
3x Fig. 3. Node numbering
sequence
for the test problems.
the Nth point, are the same as that of first point (Fig. 1). 6. The points of the boundary are to be numbered in the anti-clockwise direction [Fig. 3(a)]. If the geometry is multiply connected, the points sur-
CONCLUDING REiMARKS
For straight line segments, linear shape functions and thereby two nodes per element are sufficient. The program, however, uses quadratic shape functions uniformly. irrespective of whether the line segment is straight or curved. The two cases could have been handled separately in the program; this was not done as it would unnecessarily complicate the program and input data reading. In representing the circle, coordinates of 17 points (eight quadratic elements) have been found to be necessary [Fig. 3(a)] for computing geometric properties within 0.1% error (Table 2). On the contrary. computation using formulae would require only the coordinates of centre and radius. Therefore, in cases where the geometry is purely composed of simpler geometries such as straight lines, circles etc., the present
Table 2. Comparison of results for test problem no. I Quantity Plane lamina A xCR YCP I,, I,., I i’ Axisymmetric
x-axis
Axis of revolution: V,, St,
y - axis
J,,
J>.,
Formula
%Deviation
3.1392401 2.0000000 3.0000000 29.037384 13.341183
3.1415927 2.0000000 3 .ooooooo 29.059732 13.351769
-0.075 0.0 0.0 -0.077 -0.079
6.2809012 18.835441
6.2831853 18.849556
-0.036 -0.075
59.173282
59.217626
-0.075
4534483.3 118.39220 4243838. I
4538142.8 118.43525 4247236.2
-0.036 -0.081 - 0.080
39.448854 78.928132 1472648.9 3604395.4
39.478417 78.956835 1473926.7 3607241.7
-0.075 -0.036 -0.087 -0.079
body
Axis of revolution: V,, s J::
J,,
Program
Computing Table
3. Comparison
Quantity
geometric of results
properties for test problem
Program
863 no. 2
Formula
% Deviation
Plane lamina
A -%g )‘ce I \\ I I, I 2’
12.86075989 4.000000000 6.000000000 483.5364669 226.3212690 308.6582374 16.00000000
12.8584074 4.000000000 6.000000000 483.450600 226.282453 308.601776 16.00000000
0.018 0.0 0.0 0.018 0.017 0.018 0.0
484.8392256 603.1857895 155457137.7 144790960.0
484.7505372 603.1857895 155430997.3 144766387.9
0.018 0.0 0.017 0.017
323.2261504 402.1238597 52826940.96 121932882.1
323.1670242 402.1238597 52818808.62 121911853.4
0.018 0.0 0.015 0.017
Axisymmetric body
Axis
ofrevolution: x-axis
“,t S,,
J,, Jb,
A.xis v. ,
of revolution - y -uxis
will be inefficient. For such cases, a subroutine purely based on formulae would appear attractive. Such additional facilities could be easily added by the user. The present program, however, has been written for general purpose usage. It always handles the problems the same way, irrespective of whether the geometry is composed of simple geometric entities or not.
program
REFERENCES
1. E. Oberg and F. D. Jones, Machinery’s Handbook, 18th Edn. Industrial Press. New York (1968). 2. T. G. Hicks and S. D. Hicks (editors), Standard Handbook qf Engineering Calculations, 2nd Edn. McGraw-Hill, New York (1986).
3. H. A. Rothbart (editor), Mechanical Design and Systems Handbook, 2nd Edn. McGraw-Hill, New York (1986). 4. M. Kutz (editor), Mechanical Engineers Handbook. Wiley, New York (1986). 5. J. E. Shigley and C. R. Mischke (editors), Standard Handbook of Machine Design, McGraw-Hill, New York (1986). 6. W. C. Young, Roark’s Formulae for Stress and Strain, 6th Edn. McGraw-Hill, New York (1989). 7. E. Hinton and D. R. J. Owen, Finife Element Programming. Academic Press, London (1977). 8. 0. C. Zienkiewicz and K. Morgan, Finite Element and Approximation. Wiley, New York (1983). 9. 0. C. Zienkiewicz, Finite Element Method, 3rd Edn. Tata-McGraw-Hill, New Delhi (1979).