- Email: [email protected]
Effect of finite element mesh orientation on solution accuracy for torsional problems
Finite Elements in Analysis and Design 41 (2005) 1377 – 1383 www.elsevier.com/locate/finel
Effect of finite element mesh orientation on solution accuracy for torsional problems N. Troyania,∗ , A. Péreza , P. Baízb a Centro de Métodos Numéricos en Ingeniería, Departamento de Mecánica, Escuela de Ingeniería y Ciencias Aplicadas,
Universidad de Oriente, Venezuela b Department of Aeronautics Imperial College London, London, UK
Received 5 April 2004; accepted 31 December 2004 Available online 8 June 2005
Abstract Finite element mesh orientation has been studied in terms of the effects it may have on solution accuracy in both solid and fluid mechanics. This work studies the effect of triangular finite element mesh orientation with regards to the torsion of non-circular machine elements and structural elements in terms of maximum shearing stress computation accuracy. The membrane deflection (a Poisson problem) analogy is used to guide the automatic mesh orientation by quadrilateral diagonal swapping, starting from an initially arbitrary mesh. Two sets of results for two specific orientations are reported: mesh orientation approximately perpendicular to the membrane contour lines and mesh orientation in a direction approximately parallel to the membrane contour lines. The torsion of a square cross section member was selected as a benchmark for the analysis since it possesses a known analytic series solution that readily allows the performing of accuracy estimates. This study demonstrates that mesh orientation may have important consequences on the accuracy of the solution for the type of stress estimates treated in this work. In a more general sense, it is concluded that mesh orientation may have important consequences on the accuracy of the solution in second order Poisson type problems as well. 䉷 2005 Elsevier B.V. All rights reserved. Keywords: Finite elements; Automatic mesh orientation; Diagonal swapping; Mesh sensitivity; Torsion
∗ Corresponding author. Fax: +58 281 2697785.
E-mail address: [email protected] (N. Troyani). 0168-874X/$ - see front matter 䉷 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2004.12.009
N. Troyani et al. / Finite Elements in Analysis and Design 41 (2005) 1377 – 1383
1378
1. Introduction Finite element (FE) methods represent a powerful tool for both engineering analyses and mathematical analyses; however, they also possess certain limitations. Frequently, special solution devices can overcome these limitations. In this work, we examine the accuracy exhibited by a particular form of FE mesh manipulation when applied to the computation of maximum shearing stresses in torsional problems. To this end, a particular aspect of triangular finite element meshes is studied; namely, mesh orientation (diagonal orientation) in non-circular cross sections in order to establish the effect of this aspect on solution quality [1]. This type of mesh sensitivity analysis is of interest in both mechanical and structural engineering. Specifically, a study was performed to examine maximum shearing stress calculation accuracy for square cross section members subjected to torsion. This particular geometry was chosen since it possesses a known analytic series solution. Mesh reorientation with respect to the streamline directions was examined in fluid problems [2,3], it was also investigated in fluids flow applications with regards to variations of the directions of the velocity fields [4]. In these three works, improvements in solution accuracy are also reported. A related application for solids where orientation is considered regarding localization is reported in [5]. In the procedure that follows, the mesh orientation process is guided by the analogy that exists between the torsion of prismatic member’s problem and the deflection of a membrane problem, the latter being a Poisson type of problem, studying both the effect of aligning the FE mesh in directions approximately perpendicular to the deflected membrane contour lines (deflection =const.) and also the effect of aligning the FE mesh approximately with the direction of the deflected membrane contour lines. The latter direction is selected based on the established fact [6] that for the type of problem treated herein, the maximum stresses do occur according to these orientations. It can be shown [6] that the mathematical model of the torsion of a prismatic member is governed by the following partial differential equation (PDE) and boundary condition, j2 jx 2
+
j2 jy 2
= −2,
in R, = 0 in S
(1)
R and S represent the region of the cross section where the PDE needs to be satisfied and region boundary, respectively. (x, y) represents the Prandtl torsional function. At the same time, the mathematical model that describes the lateral deflection of a membrane attached to a flat boundary S subjected to lateral pressure p is given by the following PDE and boundary condition, j2 Z jx 2
+
j2 Z jy 2
p =− , F
in R, Z = 0 in S
(2)
z and F represent lateral deflection (generating a convex function for most of the domain R), and uniform tension per unit in-plane membrane length, respectively. R and S have the same meaning as described earlier. The two problems are mathematically analogous if Z = B , provided that B is a constant equal to p/2F . Using this analogy, an existing FORTRAN computer code [7], based on linear triangles, that calculates shearing stresses for torsional members was interfaced with a connectivity manipulating FORTRAN code, [1], to orient the mesh in an incremental procedure, as described below, based on diagonal swapping without changing the position or the numbers of the nodes from an original arbitrary FE mesh.
N. Troyani et al. / Finite Elements in Analysis and Design 41 (2005) 1377 – 1383
1379
In a related result, it was shown in [8] that linear triangles, the type of FE used in this study, exhibit better accuracy than bilinear quadrilaterals for the approximation of convex functions. In this regard, we point out that the regions of maximum shearing stress in torsion are regions where the corresponding deflected membranes are, in fact, convex. Additionally, in this reference [8] it is mentioned that mesh orientation does influence the quality of solution when using FE. The merits of using the outlined procedure are measured by comparing the results obtained herein, with those corresponding to the known analytic solution for the torsion of a square cross-sectional torsional element used as a benchmark reference.
2. Mesh orientation procedure The FE mesh orientation is based on diagonal swapping procedures as explained below and the process begins with an initial arbitrary FE mesh. Linear triangular elements were used, pairing two elements at a time to form quadrilaterals possessing two possible diagonals, Fig. 1 a and b. The deflected membrane problem is solved and corresponding contour lines are determined. At each quadrilateral centroid the level line direction is compared with the direction of the two possible diagonals. There are two possibilities: 1. Orient the diagonals according to the directions that are as close as possible to the directions perpendicular to the contour lines. 2. Orient the diagonals according to the directions that are as close as possible to the contour lines. Both cases were examined in this work. The second one is illustrated in Fig. 1. Obviously, the proposed procedure neither changes the location of nodal coordinates nor does it create more elements or nodes from the original mesh. It only changes the connectivity. Hence, it does not change the necessary computational effort at each diagonal swapping level. The procedure advances according to the following considerations and restrictions. A minimum angle to be formed in resulting triangles vertices is specified so as not to deteriorate the solution according to well-known FE results [9].
P1
P1 ‘old diagonal’ Contour line P2
P4
Contour line P2
P4
P3
P3 (a)
‘new diagonal’
(b)
Fig. 1. Swapping of diagonals to orient the mesh in the approximate directions of the contour lines. (a) Original diagonal position. (b) Swapped diagonal.
1380
N. Troyani et al. / Finite Elements in Analysis and Design 41 (2005) 1377 – 1383
The algorithm used in this work begins with an arbitrary triangulation considering all the pairs of contiguous triangles such as P1 P2 P3 and P1 P3 P4 . The swapping of the diagonals, Fig. 1, is performed if all the following conditions are satisfied: Condition 1: Quadrilateral P1 P2 P3 P4 is convex. Condition 2: The level line direction at the quadrilateral centroid strictly shows a smaller angle with the new limit P2 P4 than with the old limit P1 P3 . Condition 3: The minimal vertex angle of each of the resulting triangles P1 P2 P4 and P2 P3 P4 is strictly larger than . The pseudo-code of the procedure is analogous to the one used in [2], for a particular form of fluid flow, and is reproduced below for completeness. changes = 0 do changes = 0 for (all the triangles P1 P2 P3 ) for (for all the neighbors P1 P3 P4 ) if ((condition 1) and (condition 2) and (condition 3)) then change the diagonal of P1 P2 P3 P4 changes = changes + 1 endif endfor endfor while (changes = 0) Convergence is assured since the procedure terminates after a finite number of steps. For each set of calculations leading to a final mesh, the described procedure terminates when no more diagonal swapping takes place, implying no change in the computed shearing stresses.
3. Results and analysis A relatively low refinement mesh is used in this analysis to separate the mesh orientation aspects, as much as is reasonably possible, from the mesh refinement aspects in terms of their effect on solution accuracy. In order to illustrate the described procedure, the results corresponding to a square cross section torsional member with side size a of 0.2 m, modulus of rigidity of 77.0 GPa and subjected to an angle of twist = 5.777 × 10−03 rad/m, are presented. The known series solution to this problem is max = 60.23 MPa, obtained from the corresponding formula given in [6] and reproduced below taking the free term and infinite series first dominant term only b 8 . (3) max = a 1 − 2 sec h 2a In our case for a square cross section, b = a. Tables 1 and 2 show the resulting maximum shearing stress calculations in terms of swapping level order corresponding to a mesh of 137 nodes and 232 elements.
N. Troyani et al. / Finite Elements in Analysis and Design 41 (2005) 1377 – 1383
1381
Table 1 Maximum calculated shearing stress (in MPa) for the FE mesh orientation in a direction approximately perpendicular to the contour lines indicating swapping levels Minimum angle
5◦
10◦
15◦
20◦
25◦
30◦
40◦
50◦
Original mesh First Second Third Fourth
52.368 39.014 39.228 39.127
52.368 45.418 45.170 45.187
52.368 46.403 46.138 46.060
52.368 46.285 46.008 46.023 46.095
52.368 47.378
52.368 51.907
52.368
52.368
Table 2 Maximum shearing stress (in MPa) results for the case of aligning the mesh with the contour lines starting with resultant meshes corresponding to 5◦ and 15◦ from Table 1 Minimum angle
5◦
(From Table I for 5◦ ) Original mesh 39.127 First 52.963 Second 54.781 Third 54.759 (From Table I for 15◦ ) Original mesh 46.060 First 51.334 Second 51.657 Third 51.546 Fourth 51.488
10◦
15◦
20◦
25◦
30◦
40◦
50◦
39.127 50.822 51.178 51.220
39.127 57.665 58.333 58.317
39.127 56.412 56.611
39.127 46.540
39.127 43.265
39.127 42.129
39.127 42.101
46.060 51.697 51.626 51.589
46.060 51.608 51.758
46.060 51.432 51.475
46.060 51.605 51.593
46.060 53.529 53.565
46.060 52.301 52.315
46.060 52.506
Table 1 exhibits the results of applying the described procedure to the first of the two possibilities listed in the previous section in terms of completed levels of diagonal swapping. The initial FE calculated maximum stresses are indicated in the line labeled ‘Original Mesh’. Subsequently, the diagonal swapping procedure is carried out in order to progressively align the mesh in a direction that approximates the directions of the perpendiculars to the contour lines. The maximum stresses after the swapping procedure is completed, corresponding to each minimum angle , are the ones shown at the end of each column in boldface. It can be concluded from Table 1 that the deterioration of the quality of the results depends on the minimum angle specified. Moreover, it can be concluded that the mesh orientation affects the magnitude of the calculated shearing stresses in a significant fashion. Fig. 2a shows a typical original arbitrarily oriented FE mesh used in the study, generated by ANSYS [10]. Fig. 2b corresponds to the final FE mesh oriented in a direction approximately perpendicular to the contour lines, and Fig. 2c corresponds to the final mesh oriented in a direction approximately parallel to the contour lines. In terms of maximum shearing stress, Table 2 shows the results of applying the approximate level line mesh aligning procedure starting with the final mesh for the two cases corresponding to minimum
1382
N. Troyani et al. / Finite Elements in Analysis and Design 41 (2005) 1377 – 1383
Fig. 2. Contour lines superimposed on: (a) original arbitrarily oriented mesh; (b) final mesh oriented in a direction approximately perpendicular to the contour lines; and (c) final mesh oriented in a direction approximately parallel to the contour lines.
angles of 5◦ and 15◦ , respectively from Table 1. The resulting FE mesh and contour lines are shown in Fig. 2c. The final results for the maximum shearing stress estimates after the completion of each diagonal swapping procedure are shown in boldface for each minimum angle . Note that in some cases the results are even better than the maximum shearing stress estimate for the original arbitrarily oriented mesh, for example, 52.368 MPa in the original mesh versus 58.317 MPa in the oriented one. The observed general tendency after both procedures are applied suggests that mesh aligning in the direction perpendicular to the contour lines deteriorates the quality of the solution, whereas mesh aligning in the direction parallel to the contour lines could significantly improve maximum shearing stress calculation accuracy. The best improvement is exhibited by the results in column labeled 15◦ in Table 2 (from 5◦ ), where the original value (39.127 MPa) with respect to the analytic value (60.23 MPa) shows an error of 35.0 %, at the same time the final estimate (58.317 MPa) shows with respect to the same analytic value an error of 3.17%. It is precisely this kind of maximum shearing stress estimate difference, roughly 32.0%, which underlines the fact that mesh orientation may have significant effects on FE solution accuracy for the kind of problem treated herein.
4. Conclusions A study of mesh sensitivity in terms of mesh orientation on solution accuracy was presented showing the significance of this aspect of the FE strategy in the solution of torsional problems. The proposed study complies with basic convergence requirements since a final mesh (no more swapping possible) is reached with a finite number of swapping level procedures (no more than 4 swapping levels for our present case). Given the results exhibited in Tables 1 and 2, following conclusions can be drawn: 1. Aligning the mesh in directions parallel to the contour lines produces results that are better than the ones obtained from meshes oriented according to the directions perpendicular to the contour lines. 2. Aligning the mesh in directions parallel to the contour lines produces results that are better than the ones obtained from arbitrarily oriented meshes, for particular values of the angle . A factor to consider in actual analysis situations is that using the resulting automatically generated mesh from existing commercial codes, ones that do not consider the ideas developed in this study, could conceivably produce unfavorable mesh orientations in regions of critical stress calculations.
N. Troyani et al. / Finite Elements in Analysis and Design 41 (2005) 1377 – 1383
1383
The observations made in this analysis were also exhibited and corroborated in both elliptic and triangular torsional solid cross sections when subjected to the described mesh orientation procedures. In both cases, mesh orientation in the directions approximately parallel to the contour lines resulted in improved maximum shearing stress estimates. Finally, although this study was focused on the particular analysis of FE solutions for the maximum shearing stress determination in the torsion of solid section problems, given the existence of the analogy used in this work, it is clear that the reported observations apply to the FE solution of the general class of second order PDE problems of the Poisson type mentioned in the introductory section, as well. Acknowledgements This work was partially supported by Consejo de Investigación, Universidad de Oriente, and also supported by FUNDACITE-Anzoátegui,Venezuela. Professor Zeferino Dafonseca, Universidad del Zulia, Venezuela, kindly provided the Finite Element source code used in this work. References [1] P. Baíz, Reorientación de Elementos Finitos para la Solución de Problemas de Torsión en Elementos de Sección no Circular Basada en la Analogía de la Membrana, Tesis de Grado, Departamento de Ingeniería Mecánica, Universidad de Oriente, Venezuela, 2002. [2] T. Iliescu, A flow-aligning algorithm for convection dominated problems, Int. J. Numer. Methods Eng. 46 (7) (1999) 993–1000. [3] C.K. Chui, D. Hong, Swapping edges of arbitrary triangulations to achieve the optimal order of approximations, SIAM J. Numer. Anal. 34 (1997) 1472–1482. [4] H. Wu, I.G. Currie, Adaptive mesh refinement using a-posteriori finite element error estimation, Trans. Can. Soc. Mech. Eng. 24 (1B) (2000) 247–261. [5] P. Steinmann, K. William, Adaptive techniques for localization analysis, ASME Applied Mechanics Division, vol. 157, Adaptive, Multilevel, and Hierarchical Computational Strategies, 1992, pp. 139–164. [6] I.S. Sokolnikoff, Mathematical Theory of Elasticity, second ed., McGraw-Hill, New York, 1956. [7] Z. Da Fonseca, FORTRAN code ECONDU. FOR, Escuela de Ingeniería Mecánica, Universidad del Zulia, Venezuela, 2001. [8] E.F. D’Azevedo, Are bilinear quadrilaterals better than linear triangles?, SIAM J. Sci. Comput. 22 (1) (2000) 198–217. [9] I. Babuska, A. Aziz, On the angle condition in the finite element method, SIAM J. Numer. Anal. 13 (1976) 214–227. [10] ANSYS, user manual, version 5.4, Canonsburgh, PA, USA.
Effect of finite element mesh orientation on solution accuracy for torsional problems N. Troyania,∗ , A. Péreza , P. Baízb a Centro de Métodos Numéricos en Ingeniería, Departamento de Mecánica, Escuela de Ingeniería y Ciencias Aplicadas,
Universidad de Oriente, Venezuela b Department of Aeronautics Imperial College London, London, UK
Received 5 April 2004; accepted 31 December 2004 Available online 8 June 2005
Abstract Finite element mesh orientation has been studied in terms of the effects it may have on solution accuracy in both solid and fluid mechanics. This work studies the effect of triangular finite element mesh orientation with regards to the torsion of non-circular machine elements and structural elements in terms of maximum shearing stress computation accuracy. The membrane deflection (a Poisson problem) analogy is used to guide the automatic mesh orientation by quadrilateral diagonal swapping, starting from an initially arbitrary mesh. Two sets of results for two specific orientations are reported: mesh orientation approximately perpendicular to the membrane contour lines and mesh orientation in a direction approximately parallel to the membrane contour lines. The torsion of a square cross section member was selected as a benchmark for the analysis since it possesses a known analytic series solution that readily allows the performing of accuracy estimates. This study demonstrates that mesh orientation may have important consequences on the accuracy of the solution for the type of stress estimates treated in this work. In a more general sense, it is concluded that mesh orientation may have important consequences on the accuracy of the solution in second order Poisson type problems as well. 䉷 2005 Elsevier B.V. All rights reserved. Keywords: Finite elements; Automatic mesh orientation; Diagonal swapping; Mesh sensitivity; Torsion
∗ Corresponding author. Fax: +58 281 2697785.
E-mail address: [email protected] (N. Troyani). 0168-874X/$ - see front matter 䉷 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2004.12.009
N. Troyani et al. / Finite Elements in Analysis and Design 41 (2005) 1377 – 1383
1378
1. Introduction Finite element (FE) methods represent a powerful tool for both engineering analyses and mathematical analyses; however, they also possess certain limitations. Frequently, special solution devices can overcome these limitations. In this work, we examine the accuracy exhibited by a particular form of FE mesh manipulation when applied to the computation of maximum shearing stresses in torsional problems. To this end, a particular aspect of triangular finite element meshes is studied; namely, mesh orientation (diagonal orientation) in non-circular cross sections in order to establish the effect of this aspect on solution quality [1]. This type of mesh sensitivity analysis is of interest in both mechanical and structural engineering. Specifically, a study was performed to examine maximum shearing stress calculation accuracy for square cross section members subjected to torsion. This particular geometry was chosen since it possesses a known analytic series solution. Mesh reorientation with respect to the streamline directions was examined in fluid problems [2,3], it was also investigated in fluids flow applications with regards to variations of the directions of the velocity fields [4]. In these three works, improvements in solution accuracy are also reported. A related application for solids where orientation is considered regarding localization is reported in [5]. In the procedure that follows, the mesh orientation process is guided by the analogy that exists between the torsion of prismatic member’s problem and the deflection of a membrane problem, the latter being a Poisson type of problem, studying both the effect of aligning the FE mesh in directions approximately perpendicular to the deflected membrane contour lines (deflection =const.) and also the effect of aligning the FE mesh approximately with the direction of the deflected membrane contour lines. The latter direction is selected based on the established fact [6] that for the type of problem treated herein, the maximum stresses do occur according to these orientations. It can be shown [6] that the mathematical model of the torsion of a prismatic member is governed by the following partial differential equation (PDE) and boundary condition, j2 jx 2
+
j2 jy 2
= −2,
in R, = 0 in S
(1)
R and S represent the region of the cross section where the PDE needs to be satisfied and region boundary, respectively. (x, y) represents the Prandtl torsional function. At the same time, the mathematical model that describes the lateral deflection of a membrane attached to a flat boundary S subjected to lateral pressure p is given by the following PDE and boundary condition, j2 Z jx 2
+
j2 Z jy 2
p =− , F
in R, Z = 0 in S
(2)
z and F represent lateral deflection (generating a convex function for most of the domain R), and uniform tension per unit in-plane membrane length, respectively. R and S have the same meaning as described earlier. The two problems are mathematically analogous if Z = B , provided that B is a constant equal to p/2F . Using this analogy, an existing FORTRAN computer code [7], based on linear triangles, that calculates shearing stresses for torsional members was interfaced with a connectivity manipulating FORTRAN code, [1], to orient the mesh in an incremental procedure, as described below, based on diagonal swapping without changing the position or the numbers of the nodes from an original arbitrary FE mesh.
N. Troyani et al. / Finite Elements in Analysis and Design 41 (2005) 1377 – 1383
1379
In a related result, it was shown in [8] that linear triangles, the type of FE used in this study, exhibit better accuracy than bilinear quadrilaterals for the approximation of convex functions. In this regard, we point out that the regions of maximum shearing stress in torsion are regions where the corresponding deflected membranes are, in fact, convex. Additionally, in this reference [8] it is mentioned that mesh orientation does influence the quality of solution when using FE. The merits of using the outlined procedure are measured by comparing the results obtained herein, with those corresponding to the known analytic solution for the torsion of a square cross-sectional torsional element used as a benchmark reference.
2. Mesh orientation procedure The FE mesh orientation is based on diagonal swapping procedures as explained below and the process begins with an initial arbitrary FE mesh. Linear triangular elements were used, pairing two elements at a time to form quadrilaterals possessing two possible diagonals, Fig. 1 a and b. The deflected membrane problem is solved and corresponding contour lines are determined. At each quadrilateral centroid the level line direction is compared with the direction of the two possible diagonals. There are two possibilities: 1. Orient the diagonals according to the directions that are as close as possible to the directions perpendicular to the contour lines. 2. Orient the diagonals according to the directions that are as close as possible to the contour lines. Both cases were examined in this work. The second one is illustrated in Fig. 1. Obviously, the proposed procedure neither changes the location of nodal coordinates nor does it create more elements or nodes from the original mesh. It only changes the connectivity. Hence, it does not change the necessary computational effort at each diagonal swapping level. The procedure advances according to the following considerations and restrictions. A minimum angle to be formed in resulting triangles vertices is specified so as not to deteriorate the solution according to well-known FE results [9].
P1
P1 ‘old diagonal’ Contour line P2
P4
Contour line P2
P4
P3
P3 (a)
‘new diagonal’
(b)
Fig. 1. Swapping of diagonals to orient the mesh in the approximate directions of the contour lines. (a) Original diagonal position. (b) Swapped diagonal.
1380
N. Troyani et al. / Finite Elements in Analysis and Design 41 (2005) 1377 – 1383
The algorithm used in this work begins with an arbitrary triangulation considering all the pairs of contiguous triangles such as P1 P2 P3 and P1 P3 P4 . The swapping of the diagonals, Fig. 1, is performed if all the following conditions are satisfied: Condition 1: Quadrilateral P1 P2 P3 P4 is convex. Condition 2: The level line direction at the quadrilateral centroid strictly shows a smaller angle with the new limit P2 P4 than with the old limit P1 P3 . Condition 3: The minimal vertex angle of each of the resulting triangles P1 P2 P4 and P2 P3 P4 is strictly larger than . The pseudo-code of the procedure is analogous to the one used in [2], for a particular form of fluid flow, and is reproduced below for completeness. changes = 0 do changes = 0 for (all the triangles P1 P2 P3 ) for (for all the neighbors P1 P3 P4 ) if ((condition 1) and (condition 2) and (condition 3)) then change the diagonal of P1 P2 P3 P4 changes = changes + 1 endif endfor endfor while (changes = 0) Convergence is assured since the procedure terminates after a finite number of steps. For each set of calculations leading to a final mesh, the described procedure terminates when no more diagonal swapping takes place, implying no change in the computed shearing stresses.
3. Results and analysis A relatively low refinement mesh is used in this analysis to separate the mesh orientation aspects, as much as is reasonably possible, from the mesh refinement aspects in terms of their effect on solution accuracy. In order to illustrate the described procedure, the results corresponding to a square cross section torsional member with side size a of 0.2 m, modulus of rigidity of 77.0 GPa and subjected to an angle of twist = 5.777 × 10−03 rad/m, are presented. The known series solution to this problem is max = 60.23 MPa, obtained from the corresponding formula given in [6] and reproduced below taking the free term and infinite series first dominant term only b 8 . (3) max = a 1 − 2 sec h 2a In our case for a square cross section, b = a. Tables 1 and 2 show the resulting maximum shearing stress calculations in terms of swapping level order corresponding to a mesh of 137 nodes and 232 elements.
N. Troyani et al. / Finite Elements in Analysis and Design 41 (2005) 1377 – 1383
1381
Table 1 Maximum calculated shearing stress (in MPa) for the FE mesh orientation in a direction approximately perpendicular to the contour lines indicating swapping levels Minimum angle
5◦
10◦
15◦
20◦
25◦
30◦
40◦
50◦
Original mesh First Second Third Fourth
52.368 39.014 39.228 39.127
52.368 45.418 45.170 45.187
52.368 46.403 46.138 46.060
52.368 46.285 46.008 46.023 46.095
52.368 47.378
52.368 51.907
52.368
52.368
Table 2 Maximum shearing stress (in MPa) results for the case of aligning the mesh with the contour lines starting with resultant meshes corresponding to 5◦ and 15◦ from Table 1 Minimum angle
5◦
(From Table I for 5◦ ) Original mesh 39.127 First 52.963 Second 54.781 Third 54.759 (From Table I for 15◦ ) Original mesh 46.060 First 51.334 Second 51.657 Third 51.546 Fourth 51.488
10◦
15◦
20◦
25◦
30◦
40◦
50◦
39.127 50.822 51.178 51.220
39.127 57.665 58.333 58.317
39.127 56.412 56.611
39.127 46.540
39.127 43.265
39.127 42.129
39.127 42.101
46.060 51.697 51.626 51.589
46.060 51.608 51.758
46.060 51.432 51.475
46.060 51.605 51.593
46.060 53.529 53.565
46.060 52.301 52.315
46.060 52.506
Table 1 exhibits the results of applying the described procedure to the first of the two possibilities listed in the previous section in terms of completed levels of diagonal swapping. The initial FE calculated maximum stresses are indicated in the line labeled ‘Original Mesh’. Subsequently, the diagonal swapping procedure is carried out in order to progressively align the mesh in a direction that approximates the directions of the perpendiculars to the contour lines. The maximum stresses after the swapping procedure is completed, corresponding to each minimum angle , are the ones shown at the end of each column in boldface. It can be concluded from Table 1 that the deterioration of the quality of the results depends on the minimum angle specified. Moreover, it can be concluded that the mesh orientation affects the magnitude of the calculated shearing stresses in a significant fashion. Fig. 2a shows a typical original arbitrarily oriented FE mesh used in the study, generated by ANSYS [10]. Fig. 2b corresponds to the final FE mesh oriented in a direction approximately perpendicular to the contour lines, and Fig. 2c corresponds to the final mesh oriented in a direction approximately parallel to the contour lines. In terms of maximum shearing stress, Table 2 shows the results of applying the approximate level line mesh aligning procedure starting with the final mesh for the two cases corresponding to minimum
1382
N. Troyani et al. / Finite Elements in Analysis and Design 41 (2005) 1377 – 1383
Fig. 2. Contour lines superimposed on: (a) original arbitrarily oriented mesh; (b) final mesh oriented in a direction approximately perpendicular to the contour lines; and (c) final mesh oriented in a direction approximately parallel to the contour lines.
angles of 5◦ and 15◦ , respectively from Table 1. The resulting FE mesh and contour lines are shown in Fig. 2c. The final results for the maximum shearing stress estimates after the completion of each diagonal swapping procedure are shown in boldface for each minimum angle . Note that in some cases the results are even better than the maximum shearing stress estimate for the original arbitrarily oriented mesh, for example, 52.368 MPa in the original mesh versus 58.317 MPa in the oriented one. The observed general tendency after both procedures are applied suggests that mesh aligning in the direction perpendicular to the contour lines deteriorates the quality of the solution, whereas mesh aligning in the direction parallel to the contour lines could significantly improve maximum shearing stress calculation accuracy. The best improvement is exhibited by the results in column labeled 15◦ in Table 2 (from 5◦ ), where the original value (39.127 MPa) with respect to the analytic value (60.23 MPa) shows an error of 35.0 %, at the same time the final estimate (58.317 MPa) shows with respect to the same analytic value an error of 3.17%. It is precisely this kind of maximum shearing stress estimate difference, roughly 32.0%, which underlines the fact that mesh orientation may have significant effects on FE solution accuracy for the kind of problem treated herein.
4. Conclusions A study of mesh sensitivity in terms of mesh orientation on solution accuracy was presented showing the significance of this aspect of the FE strategy in the solution of torsional problems. The proposed study complies with basic convergence requirements since a final mesh (no more swapping possible) is reached with a finite number of swapping level procedures (no more than 4 swapping levels for our present case). Given the results exhibited in Tables 1 and 2, following conclusions can be drawn: 1. Aligning the mesh in directions parallel to the contour lines produces results that are better than the ones obtained from meshes oriented according to the directions perpendicular to the contour lines. 2. Aligning the mesh in directions parallel to the contour lines produces results that are better than the ones obtained from arbitrarily oriented meshes, for particular values of the angle . A factor to consider in actual analysis situations is that using the resulting automatically generated mesh from existing commercial codes, ones that do not consider the ideas developed in this study, could conceivably produce unfavorable mesh orientations in regions of critical stress calculations.
N. Troyani et al. / Finite Elements in Analysis and Design 41 (2005) 1377 – 1383
1383
The observations made in this analysis were also exhibited and corroborated in both elliptic and triangular torsional solid cross sections when subjected to the described mesh orientation procedures. In both cases, mesh orientation in the directions approximately parallel to the contour lines resulted in improved maximum shearing stress estimates. Finally, although this study was focused on the particular analysis of FE solutions for the maximum shearing stress determination in the torsion of solid section problems, given the existence of the analogy used in this work, it is clear that the reported observations apply to the FE solution of the general class of second order PDE problems of the Poisson type mentioned in the introductory section, as well. Acknowledgements This work was partially supported by Consejo de Investigación, Universidad de Oriente, and also supported by FUNDACITE-Anzoátegui,Venezuela. Professor Zeferino Dafonseca, Universidad del Zulia, Venezuela, kindly provided the Finite Element source code used in this work. References [1] P. Baíz, Reorientación de Elementos Finitos para la Solución de Problemas de Torsión en Elementos de Sección no Circular Basada en la Analogía de la Membrana, Tesis de Grado, Departamento de Ingeniería Mecánica, Universidad de Oriente, Venezuela, 2002. [2] T. Iliescu, A flow-aligning algorithm for convection dominated problems, Int. J. Numer. Methods Eng. 46 (7) (1999) 993–1000. [3] C.K. Chui, D. Hong, Swapping edges of arbitrary triangulations to achieve the optimal order of approximations, SIAM J. Numer. Anal. 34 (1997) 1472–1482. [4] H. Wu, I.G. Currie, Adaptive mesh refinement using a-posteriori finite element error estimation, Trans. Can. Soc. Mech. Eng. 24 (1B) (2000) 247–261. [5] P. Steinmann, K. William, Adaptive techniques for localization analysis, ASME Applied Mechanics Division, vol. 157, Adaptive, Multilevel, and Hierarchical Computational Strategies, 1992, pp. 139–164. [6] I.S. Sokolnikoff, Mathematical Theory of Elasticity, second ed., McGraw-Hill, New York, 1956. [7] Z. Da Fonseca, FORTRAN code ECONDU. FOR, Escuela de Ingeniería Mecánica, Universidad del Zulia, Venezuela, 2001. [8] E.F. D’Azevedo, Are bilinear quadrilaterals better than linear triangles?, SIAM J. Sci. Comput. 22 (1) (2000) 198–217. [9] I. Babuska, A. Aziz, On the angle condition in the finite element method, SIAM J. Numer. Anal. 13 (1976) 214–227. [10] ANSYS, user manual, version 5.4, Canonsburgh, PA, USA.