Three-dimensional mesh generation using principles of finite element method

Advances in Engineering Software 31 (2000) 25–34 www.elsevier.com/locate/advengsoft

Three-dimensional mesh generation using principles of finite element method M. Karumbu Nathan, S. Prasanna, G. Muthuveerappan* Machine Design Section, Department of Mechanical Engineering, Indian Institute of Technology, Madras, Chennai 36, India Received 1 December 1998; received in revised form 7 April 1999; accepted 11 May 1999

Abstract The present work deals with the development of a three-dimensional mesh generation algorithm using the principles of FEM with special emphasis on the computational efficiency and the memory requirement. The algorithm makes use of a basic mesh that defines the total number of elements and nodes. Wavefront technique is used to renumber the nodes in order to reduce the bandwidth. By elastic distortion of the basic mesh, it is redefined to map onto actual geometry to be discretized. Later a finer distribution of mesh is done in the zones of interest to suit the nature of the problem. The same Finite Element code meant for stress analysis is adopted with necessary modifications. The algorithm has been extended to three-dimensional geometries. The current methodology is used to discretize a straight bevel gear and an hourglass worm to study their stress patterns. q 2000 Elsevier Science Ltd. All rights reserved. Keywords: Automatic mesh generation; Basic mesh; FEM; Elastic distortion; Influence coefficient; Wavefront technique

1. Introduction FEM has come to stand as a very powerful tool in static and dynamic analysis of many engineering problems. The applications of FEM to certain problems like stress analyses in gears become difficult due to complicated geometry of the tooth gearing. Moreover, the stresses evaluated are very much influenced by the type of mesh used to discretize the gear geometry. Making an appropriate element division and determination of the nodal co-ordinates, for the purpose of stress analysis, become very cumbersome. Apart from these factors, one has to come across zones of stress concentration in the gear geometry where a fine distribution of elements is necessary. Hence a mesh generation algorithm should include all features to overcome these problems. The generation of a Finite Element mesh is a geometrically controlled process. Various schemes are available to perform the task of mesh generation. A method involving two stages to generate the mesh automatically was proposed by Buell et al. [1]. Here the intermediate nodes were first interpolated from the boundary nodes and then connected to form the elements. Most of the literature in the area of automatic mesh generation was focused on the generation of triangular and tetrahedral elements [2–4]. The specified * Corresponding author. Tel.: 1 91-044-235-1365; fax: 1 91-044-2350509. E-mail address: [email protected] (G. Muthuveerappan)

region is discretized using triangular elements. An efficient algorithm [5] for Delaunay triangulation was presented in 1996 and many acceleration procedures were implemented to speed up the process. However, there involved a number of critical issues in the development of a Delaunay-based automatic mesh generator. In cases of stress analysis, users have a strong bias against the use of triangular and tetrahedral elements. The constant strain triangular elements behave poorly in bending as its strain field contains only constants. The second order tetrahedral elements provide a cost effective means for performing stress analysis. In many two- and three-dimensional analyses, use of quadrilateral and hexahedral elements are preferred over other types of elements. In 1991, a new approach called Paving Theory [6], was proposed wherein the region was divided using quadrilateral elements. There was another class of algorithm that made use of parametric mapping [7]. The advantage of the mapping method over triangulation scheme is that, it does not require extensive geometric checking for connecting the nodes in elements. The element size and density could be varied by suitably varying the successive parameters at the boundary. In 1988, the status of methods [8] for automatic generation and control of Finite Element meshes were presented with much emphasis on three-dimensional mesh generation. An updated scenario in the development of automatic mesh generation for three-dimensional domains was presented in 1996 by the same author [9]. In his work, various means of

0965-9978/00/$ - see front matter q 2000 Elsevier Science Ltd. All rights reserved. PII: S0965-997 8(99)00024-1

26

M. Karumbu Nathan et al. / Advances in Engineering Software 31 (2000) 25–34

(b)

(a)

Fig. 1. Example to show the selection of a basic mesh: (a) circles; (b) 8-sided polygon.

improving the Finite Element mesh to attain the desired level of accuracy was provided, namely, by either subdividing the selected elements into smaller elements of the same type (h-refinement) or by increasing the order of polynomial of selected elements (p-refinement) or by relocating the position of nodes within a given Finite Element mesh topology (r-refinement). A very simple algorithm to generate the mesh automatically is proposed in the present work using Finite Element Method. This algorithm is similar to the process of r-refinement where the nodes are relocated at each level of refinement. As Finite Element procedures and solution techniques are very familiar, a method that incorporates all the features of the conventional finite element approach is used to generate the mesh. The advantage of the present method is that, it is easy to program. The same finite element source code can be used to generate the mesh and subsequently to analyse the structure. Special emphasis is laid to improve the efficiency of the process and to reduce the computational labour. An algorithm to discretize a spur gear model in a two-dimensional plane using triangular elements was

Fig. 2. Basic mesh.

proposed by Suzuki et al. [10], however, due importance was not given to the efficiency of the program in terms of computational labour. The current work is an extension of this algorithm with added improvements to handle threedimensional geometries. In the first step of the generation of a mesh, the domain of the model is defined. Based on the domain a valid basic mesh is selected with optimum number of nodes and elements. The nodal co-ordinates of the intermediate nodes are determined by using the same Finite Element routine, with the only change being the incorporation of a new elemental stiffness matrix. Keeping in view the size of problem involved, an algorithm is developed for renumbering the nodes based on wavefront technique, thus reducing core and computational labour involved. In this work the mesh generation scheme is used to discretize an hourglass worm gear and a straight bevel gear. In the zones of interest, a fine distribution of mesh is done. To validate the mesh generated, the stress patterns in the gears are studied by applying random loads.

2. Concept of basic mesh In the mesh generation process, a basic mesh consisting of a set of regular meshes should be created first. This becomes the basis for analysis of the model. The peripheral surface geometry of a model should always correspond to a good basic mesh that is always a domain comprising of straight edges at the boundary. The optimum number of nodes defining the shape of the basic mesh should be derived such that it is sufficient to represent geometric irregularity of boundary of the model. For example, the eightsided polygon as shown in Fig. 1(b) can be assumed as the shape of basic mesh to represent the annular region between the circles shown in Fig. 1(a). Each side of the polygon is further divided using number of nodes. However, a better representation of the model is achieved if number of sides of the polygon is increased. The basic mesh can be considered to be a rubber net stretched over a surface whose geometry is well defined. The line segment connecting each node to its adjacent node is assumed to be a linear elastic member and this is referred as influence coefficient. The basic mesh becomes distorted in a manner such that boundary nodes are made to coincide with the exact boundary of the model under analysis. The reorientation of interior nodes of the basic mesh gives co-ordinates of interior nodes of the actual model. In this process, initially it is assumed to have equal influence coefficient for all line segments connecting the nodes. To get a finer distribution of mesh, the influence coefficients can be varied. Special emphasis is laid on the selection of a proper basic mesh as the discretized geometry is entirely dependent on the basic mesh. A basic mesh, shown in Fig. 2, consists of a set of regular quadrilateral elements. This basic mesh is used to discretize an hourglass worm shown in Fig. 3 in its axial plane. Each

M. Karumbu Nathan et al. / Advances in Engineering Software 31 (2000) 25–34

m

X1

27

n

X2

X3

Fig. 5. Example of one-dimensional model to demonstrate the nodal co-ordinate calculations based on m:n.

manner such that its boundaries coincide with the periphery of the exact tooth profile. 3. Principle of mesh generation

Fig. 3. A discretized hourglass worm with equal influence coefficients for all the line segments.

intermediate node in Fig. 3 is placed at the centroid of nodes that are connected to it, i.e. every line segment connecting the nodes is given equal influence coefficient. The shape of the mesh attained by varying the influence coefficients of line segments near the fillet region is shown in Fig. 4. This procedure is similar to increasing the elasticity of line segments near the fillet region as a result of which the nodes near the fillet are reoriented to get a finer mesh distribution. The nodes on the boundary are referred as fixed nodes. The interior nodes to be determined for this model are referred as free nodes. This basic mesh is stretched in a

The co-ordinates of all interior nodes have to be determined exactly, based on co-ordinates of the boundary nodes, ensuring that the resultant mesh is well balanced for FEM calculations. Considering the case of a one-dimensional element shown in Fig. 5, X1, X2 represent the co-ordinates of ends of a line segment and X3 is the co-ordinate of the node which divides the segment in the ratio m:n. The co-ordinate of X3 is defined by     m n X3 ˆ X 1 X : m1n 2 m1n 1 The co-ordinate of X3 solely depends on the ratio m:n. By reducing m, the value of X3 is brought closer to X1 and vice versa. By defining K1 and K2 such that m:n is equal to K2:K1, the value of X3 is obtained by the following equation:     K2 K1 X3 ˆ X 1 X K1 1 K2 2 K1 1 K2 1 i.e. 2K1 X1 1 …K1 1 K2 †X3 2 K2 X2 ˆ 0 where Kis are the influence coefficient of line segments. Extending this concept further, it can be applied to another one-dimensional problem shown in Figs. 6 and 7 with two boundary nodes (fixed nodes) and three undefined interior nodes (free nodes). Their respective equations can be written as: …K1 1 K2 †X2 2 K2 X3 ˆ K1 X1 2K2 X2 1 …K2 1 K3 †X3 2 K3 X4 ˆ 0 2K3 X3 1 …K3 1 K4 †X4 ˆ K4 X5 :

K1 X1 Fig. 4. An example of discretization with varying influence coefficients at the zones of interest (fillet).

K2

X3

X2

Fig. 6. Example of one-dimensional model to demonstrate the nodal co-ordinate calculations using influence coefficients.

28

M. Karumbu Nathan et al. / Advances in Engineering Software 31 (2000) 25–34

Fig. 7. Example of one-dimensional model to demonstrate the nodal co-ordinate calculations.

These equations can be put in a matrix form with X2, X3 and X4 as the unknowns. To get a variation in the distribution of the intermediate nodes the influence coefficients K1, K2, K3, K4 can be varied as per the requirement: 9 2 38 9 8 2K2 0 K1 1 K2 X2 > K1 X1 > > > > > > > = 6 7< = < 6 2K2 7 K2 1 K3 2K3 5 X3 ˆ 0 : 4 > > > > > > : > : ; > ; 0 2K3 K3 1 K4 X4 K4 X5 The LHS of the above matrix is similar to the assembled stiffness matrix in FEM with the Kis representing the influence coefficient of line segments. The RHS of the above matrix is similar to the force vector. Solving the above equations the intermediate co-ordinates can be determined. This procedure can easily be extended to two-dimensional as well as three-dimensional problems. Since this algorithm is similar to the assembly of stiffness matrices and computation of resultant deflections, as used in conventional FEM, the same code can be used with a little modification for the purpose of mesh generation. In this algorithm the stiffness matrix of each line segment is given by " # K 2K 2K

K

where K is the influence coefficient of the line segment. Based on the connectivity details, the stiffness matrix of each line segment is assembled to get the global stiffness matrix. The force vector entirely depends on the product of fixed nodes and influence coefficient of line segments

Fig. 8. Node numbering before applying wavefront technique.

connecting the corresponding fixed nodes. With this small change incorporated in the conventional FEM, any geometrical model can be discretized. The Kis are initially assigned a value of unity, i.e. each node is placed at the centroid of its surrounding nodes that are connected to it. As the refinement is required in certain regions, the values of Ki in that region alone are altered. It may appear that the variation of influence coefficient of line segments in one region can distort the element elsewhere. However, it could be observed from the assembled stiffness matrix that, the variation of influence coefficients, Kis will be more predominantly felt in the zone where they have been altered. The discrepancies in the element irregularity will be projected when an error analysis is performed to check for the aspect ratio, interior angle, etc. Once the particular portion of the mesh requiring refinement is determined, the values of Ki in that region are changed accordingly.

4. Wave front technique Keeping in view the size of mesh involved in threedimensional problems, a suitable numbering scheme should be adopted so as to attain optimum bandwidth. Various numbering schemes are available for optimising the bandwidth. Popular among them are Collins bandwidth reduction scheme [11] and Cuthill–McKee scheme [12]. The disadvantage of Collins algorithm is that it cannot be applied to 8-noded hexahedral elements. Likewise in the Cuthill–McKee scheme the efficiency of the algorithm depends entirely on the selection of the starting node for the numbering process. A code is developed in the present work to renumber the nodes based on the wavefront technique [13]. In this technique importance is given to element numbering rather than node numbering. This scheme has an advantage over the other schemes in terms of computational time involved to attain satisfactory nodal numbering. Assuming a wave to sweep through the surface to be discretized, the elements are numbered sequentially along the direction of propagation of the wave. The direction of propagation should be such that the wave passes through minimum number of elements. It is found that when the wave propagates through more than two elements at a time, zigzag order of element numbering is found to be more efficient. The first step in this numbering scheme is to assign the element numbers, based on the principle of wave propagation explained previously. The element numbers remain constant in this scheme while the nodes are renumbered. Initial nodal numbering is done in a random order but when it is proceeded from one element to the other, the node that is eliminated first is replaced as the first node. By adopting this procedure, terms involved in the stiffness matrix are brought closer to the principal diagonal. As a

M. Karumbu Nathan et al. / Advances in Engineering Software 31 (2000) 25–34

29

Table 1 Details of node numbers after the application of wavefront technique with reference to Figs. 8 and 9 (core required before renumbering ˆ 141 (in terms of the space required to store the data), core required after renumbering ˆ 104)

Fig. 9. Node numbering after applying wavefront technique.

typical example, in the structure shown in Figs. 8 and 9, when it is discretized using triangular elements, the numbers (1), (2), etc. represent the element numbers while 1, 2, etc. refer to the node numbers. The nodes are initially numbered in a random order. Table 1 shows the order in which the nodes are being eliminated, and the corresponding new node numbers.

5. Concept of three-dimensional mesh generation The concept of generation of intermediate nodes in a onedimensional line segment can be extended to two as well as three-dimensional structures. Assemblage of stiffness matrix is based on the fact that each element is a line segment. The generation of connectivity of a large structure as an assemblage of individual line segments is more cumbersome and hence a new concept is proposed to overcome this problem. The nodal connectivity generated for stress analysis can be used for the purpose of mesh generation to avoid the labour of giving the connectivity as individual line segments. From the hexahedral element shown in Fig. 10, it can be inferred that each element is an assemblage of twelve line segments whose stiffness matrix (8 × 8) is given by 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4

Element

End1

End2

End3

Node getting eliminated

Corresponding new node number

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

10 10 5 5 6 6 11 11 12 12 13 13 15 15 20 20 25 25 24 24 23 23 22 22

4 5 1 2 3 2 6 5 6 7 8 7 13 12 15 14 20 19 19 18 18 17 17 16

9 4 4 1 2 5 5 10 11 6 7 12 12 14 14 19 19 24 18 23 17 22 16 21

9 – 4 1 3 2 – 5,10 11 6 8 7 13 12 15 14 20 25 19 24 18 23 17 22,16,21

1 – 2 3 4 5 – 6,7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23,24,25

The above stiffness matrix is the basic elemental stiffness matrix. Similarly, if tetrahedral elements are used to discretize a structure, the element can be considered as an assemblage of six line segments as shown in Fig. 11 and the corresponding elemental stiffness matrix is given by 2 6 6 6 6 6 6 4

K4 1 K5 1 K6

2K4

2K5

2K6

2K4

K1 1 K2 1 K4

2K1

2K5

2K5

2K1

K1 1 K3 1 K5

2K3

2K6

2K5

2K3

K3 1 K5 1 K6

3 7 7 7 7: 7 7 5

The entire structure can either be expressed as an assemblage of (8 × 8) elemental stiffness matrix if hexahedral

K1 1 K4 1 K8

2K1

0

2K4

2K8

0

0

0

2K1

K1 1 K2 1 K5

2K2

0

0

2K5

0

0

0

2K2

K2 1 K3 1 K6

2K3

0

0

2K6

0

2K4

0

2K3

K3 1 K4 1 K7

0

0

0

2K7

2K8

0

0

0

K8 1 K9 1 K12

2K9

0

2K12

0

2K5

0

0

2K9

K5 1 K9 1 K10

2K10

0

0

0

2K6

0

0

2K10

K6 1 K10 1 K11

2K11

0

0

0

2K7

2K12

0

2K11

K7 1 K11 1 K12

3 7 7 7 7 7 7 7 7 7 7 7: 7 7 7 7 7 7 7 7 5

30

M. Karumbu Nathan et al. / Advances in Engineering Software 31 (2000) 25–34

Fig. 10. 8-noded hexahedral element.

element or as an assemblage of (4 × 4) elemental stiffness matrix if tetrahedral element is used. For any other element, the elemental stiffness matrix can be derived using the same argument. It is to be noted that when an influence coefficient is assigned to one line segment, it has to be distributed equally to all the adjacent hexahedral elements. The source code developed uses skyline approach to compute the nodal co-ordinates. Computational labour in terms of the number of multiplication’s involved and core required for mesh generation are calculated using the above code. Three different kinds of meshes generated from the same basic mesh (Fig. 12) with equal influence coefficient of line segments are shown in Figs.13–15. The examples are illu-

Fig. 11. 4-noded tetrahedral element.

Fig. 12. Basic mesh of a three dimensional solid.

Fig. 13. Example of the generated mesh showing the elastic distortion for a three-dimensional solid. (All nodes on the faces constrained.)

M. Karumbu Nathan et al. / Advances in Engineering Software 31 (2000) 25–34

31

of the cube are alone constrained and the resulting mesh is shown in Fig. 15.

6. Outline of the mesh generation algorithm

Fig. 14. Example of the generated mesh showing the elastic distortion for a three-dimensional solid. (All nodes on edges constrained.)

strated to show the effect of proper constraints on the resultant mesh. If all the nodes on six faces are constrained, the resulting mesh obtained would be as shown in Fig. 13. The resulting mesh shown in Fig. 14 is obtained by constraining the 12 edges of the cube. It is observed that the intermediate nodes on the faces are drawn towards the centre. The reason for this behaviour is that the nodes on the faces have influence of the interior nodes apart from the fixed nodes. This behaviour is very predominant when the eight corner nodes

The steps involved in the mesh generation algorithm is explained in the flowchart shown in Fig. 16. The first step in the generation of mesh involves the development of a basic mesh based on the contour of the model to be discretized. The basic mesh, is then divided into a set of regular meshes using optimum number of elements and nodes. The set of divisions constituting the basic mesh gives the topological relationship between the nodes. The shape of the basic mesh should be based on the irregularity of boundary of the model. Sufficient number of nodes should be selected on boundary of the basic mesh so that the resultant basic mesh represents the boundary of the actual model. The nodal connectivity developed for the purpose of stress computation is used as an input for the purpose of mesh generation. The influence coefficients, Kis, of all line segments are initially assigned a value of unity. The nodes are then renumbered based on wavefront technique to optimize the bandwidth. The initial value of Kis is used to determine the elemental stiffness matrix. The elemental stiffness matrices are then assembled based on the connectivity details and the nodal co-ordinates are then determined. An error analysis is finally performed to check for the validity of the element. The influence coefficients of line segments are varied in the zones of interest where further refinement of the mesh is necessary. Refinement of the mesh could also be performed by varying the location of fixed nodes on the external contour of the model.

7. Application of the algorithm

Fig. 15. Example of the generated mesh showing the elastic distortion for a three-dimensional solid. (All corner nodes constrained.)

Gears have complicated geometry with stress concentration usually predominant in the fillet region. Therefore a very fine distribution of mesh is required in the fillet region. To validate the scheme, two types of gear teeth have been discretized. The geometries of gears are based on the specifications provided by Dudley [14]. In the first case, the three-dimensional mesh generation algorithm is used to discretize a straight bevel gear tooth with fine distribution of elements around the fillet region. In view of the magnitude of the problem for stress analysis, a minimum number of elements, without loss of accuracy of the results, is selected. A total of 510 elements and 762 nodes are used to discretize the gear tooth. The outline of the gear tooth before discretization and the corresponding basic mesh is given in Figs. 18 and 17, respectively. The mesh of the tooth after discretization is shown in Fig. 19. The hidden traces of the back portion of the tooth are removed and this is shown in Fig. 20. The time

32

M. Karumbu Nathan et al. / Advances in Engineering Software 31 (2000) 25–34

Fig. 16. Various steps involved in mesh generation.

taken to discretize the bevel gear tooth in IBM RS-6000 is 6 s. In the second case a portion of the hourglass worm thread ( ^ 108 on either sides of the axial plane), shown in Fig. 22 , is discretized in the same manner as explained earlier using a similar basic mesh (Fig. 21). The basic meshes of the bevel gear and the hourglass worm are generated to take care of irregularity in the boundary of the actual models. For better clarity of the worm an additional view (Fig. 23) is shown after hiding the traces of back portions of the worm. In

actual analysis a total of 660 elements and 957 nodes are taken into consideration ( ^ 508 on either side of the axial plane). To prove the validity of the discretized model of worm and the bevel gears, a check for convergence of stresses is performed by loading the tooth in each case. The mesh is finalised after convergence of stresses is attained. A good discretized model of the tooth is attained after four iterations by varying the influence coefficients of the line segments.

M. Karumbu Nathan et al. / Advances in Engineering Software 31 (2000) 25–34

33

Fig. 19. Automatically generated mesh of a straight bevel gear. Fig. 17. Basic mesh for a straight bevel gear.

8. Conclusion An improved algorithm over the one proposed by Suzuki [10], is presented in this work involving the features of Finite element Method. The salient features of the algorithm are as follows:

reducing the core required for storing data by this technique. 5. The main feature of the program is that the same finite element analysis code can be used for mesh generation with minor modifications in the stiffness matrix. 6. Based on the desired accuracy of results, the mesh is made finer. This procedure can be automated thereby making this algorithm adaptive.

1. Easy generation of input data. 2. Refinement of mesh can be automatically done by changing influence coefficients of line segments. 3. Extension of the algorithm to three-dimensional objects is easier, as the same nodal connectivity defined for stress analysis purpose using FEM can be used. 4. Computational time has been minimized by effectively

The procedure of mesh generation scheme presented in this work has certain difficulties as it does not introduce additional elements or nodes into the system. Hence the accuracy of results is limited to a certain extent which is governed by number of elements used for dicretization.

Fig. 18. External contour of a straight bevel gear tooth.

Fig. 20. Automatically generated mesh of a straight bevel gear with hidden traces removed.

34

M. Karumbu Nathan et al. / Advances in Engineering Software 31 (2000) 25–34

Fig. 23. Automatically generated mesh of an hourglass worm with hidden traces removed. Fig. 21. Basic mesh for the hourglass worm.

Apart from that automatic generation of basic mesh for complicated geometries is difficult. The authors are currently working in these areas to overcome the above mentioned problems.

Acknowledgements The authors are very grateful to Prof V. Ramamurti for the valuable suggestions in executing this work.

Fig. 22. Automatically generated mesh of an hourglass worm.

References [1] Buell WR, Bush BA. Mesh generation—a survey, transactions of ASME. Journal of Engineering for Industry 1973;:332–338. [2] Bowyer A. Computing Dirichlet tesselations. The Computer Journal 1981;24(21):162–166. [3] Boender E. Reliable Delaunay-based mesh generation and mesh improvement. Communications in Numerical methods in Engineering 1994;10:773–783. [4] Joe B, Simpson RB. Triangular meshes for regions of complicated shape. International Journal for Numerical Methods in Engineering 1986;23:751–778. [5] Borouchaki H, George PL. Optimal Delaunay point insertion. International Journal for Numerical Methods in Engineering 1996;39:3407–3437. [6] Blacker D, Stephenson MB. Paving: a new approach to automated quadrilateral mesh generation. International Journal for Numerical Methods in Engineering 1991;32:811–847. [7] Crawford H, Waggenspack WN, Anderson DC. Composite mapping for planar mesh generation. International Journal for Numerical Methods in Engineering 1987;24:2241–2252. [8] Shephard MS. Approaches to the automatic generation and control of finite element meshes. Applied Mechanics Reviews 1988;41(4):169– 185. [9] Shephard MS. Update to: approaches to the automatic generation and control of finite element meshes. Applied Mechanics Reviews, vol. 49, no 10, part 2, October 1996. pp S5–S14. [10] Suzuki T, Aida T, Kubo A, Chong TH, Fujio H. Tooth fillet stresses of gear with thin rim. Bulletin of the Japanese Society of Mechanical Engineers 1982;:1022–1029. [11] Collins RJ. Bandwidth reduction by automatic renumbering. International Journal for Numerical Methods in Engineering 1973;6:345– 356. [12] Cook RD. Concepts and applications of finite element analysis, New York: Wiley, 1981. [13] Ramamurti V. Computer aided mechanical design and analysis, 4. New York: McGraw-Hill, 1998. [14] Dudley DW. Practical gear design, New York: McGraw-Hill, 1954.