Comparison of calculated stress; figures in parenthesis indicate a rough estimate of computational time
1. Introduction
In the last decades, finite element (FE) method has proven to be an efficient tool for the numerical analysis of two- or three-dimensional structures of whatever complexity, in mechanical, thermal or other physical problems. It is widely recognised that computational cost (as time and computer capability) increases greatly with structure complexity, being larger with three-dimensional analyses than with two-dimensional ones.
It is therefore desirable to devise simplified approaches that may provide a reduction in overall computational effort. An example of considerable importance is the study of bodies of revolution (axisymmetric structures) under axisymmetric loading, where a three-dimensional problem is solved by a two-dimensional analysis. Examples are vessels under internal pressure, rotating disks, foundation piles (Zienkiewicz & Taylor, 2000).
Apart from axisymmetric problems solved by a plane model, a full three-dimensional analysis is needed, in principle, whenever the structure is axisymmetric but the load is not. In such situations, often encountered in many engineering applications, it is desirable to search for simplified approaches, which may still replace (and thus avoid the computation effort needed by) the use of full three-dimensional simulations.
A particular sub-class of problems is encountered when the load applied to axisymmetric structures is exactly antisymmetric; an example is a shaft under a torsion load, which, as it will be shown, can be solved by a plane FE approach, which greatly simplifies the analysis.
Another example is represented by semi-analytical methods, which have been developed more than fifty years ago for FE analysis of axisymmetric structures loaded non-axisymmetrically (Wilson, 1965). Such methods use a Fourier series expansion to reduce a three-dimensional problem to a two-dimensional harmonic model and to compute the solution as superposition of results of every harmonic component analysis. At present, this approach is still not well established and it is rarely used in mechanical design. Even commercial FE codes including harmonic elements have found limited application, due to practical difficulties related to Fourier series conversion of external loads. Only few applications of semi-analytical methods to engineering practical cases have been reported in literature, see (Genta & Tonoli, 1996; Lai & Booker, 1991; Kim et al., 1994; Pedersen & Laursen, 1982; Taiebait & Carter, 2001; Thomas et al., 1983; Zienkiewicz & Taylor, 2000).
The goal of this work is twofold. First, it aims to provide a theoretical background on the use of semi-analytical FE approach in numerical analysis of axisymmetric structures loaded non-axialsymmetrically. In particular, two original results are developed: a plane "axi-antisymmetric" FE model for solving axisymmetric components loaded in torsion, a semi-analytical approach for the analysis of plane axisymmetric bodies under non-axisymmetric thermal loadings. The second intent of the work is to clarify some practical aspects in the application of semi-analytical method to engineering problems. Two illustrative examples are then discussed: an axisymmetric component with shoulder fillet under axial, bending and torsion mechanical loading, followed by the more general case of a shaft; a simplified numerical approach for estimating the transient temperature field in a rotating cylinder under thermal loadings. The presented examples confirm how the proposed approach gives a high accuracy in results, with also a significant reduction in computational time, compared to classical FE analyses.
2. Theoretical background
A three-dimensional structure or solid is defined as "axisymmetric" if its geometry, material properties and boundary conditions are independent of an azimuth coordinate
Depending on the configuration of external loads, different types of analysis can be identified. For example, if also external loads are themselves axisymmetric with respect to same
Another situation occurs for an axisymmetric structure under an "axi-antisymmetric" loading (this term will be clarified later on), i.e. a loading which is axial anti-symmetric with respect to
A third situation, of great practical interest, occurs when a structure is axisymmetric but the loading is not, so that the analysis becomes really three-dimensional. A great simplification can be achieved by using a so-called semi-analytical approach, which adopts a harmonic model based on Fourier series method.
The next sections will address the main aspects of FE theory for the case of two- and three-dimensional axisymmetric structures under different types of both mechanical and thermal loadings. For a mechanical analysis, the theory for axisymmetric, axi-antisymmetric and non-axisymmetric mechanical loads will be developed. For a thermal analysis, a one-dimensional finite element for dealing with axisymmetric plane body under non-axisymmetric thermal loadings will be next developed; the steady-state and transient case will be considered.
3. Mechanical analysis
This section is concerned with FE theory for mechanical analysis of axisymmetric structures under, respectively, axisymmetric, axi-antisymmetric and non-axisymmetric loads.
3.1. Plane axisymmetric finite element
Although the theory of axisymmetric mechanical analysis is well known (Cook et al., 1989; Wilson, 1965; Zienkiewicz & Taylor, 2000), it will be shortly reviewed to introduce several equations that will be used in the following sections. In the case of axisymmetric structures loaded by axially symmetric loads, by symmetry, the two displacement components
where
As already mentioned, four non-zero strain components have to be considered in an axisymmetric deformation, see Fig. 2; the strain vector in polar coordinates thus is:
where submatrix [
where
The element stiffness matrix is calculated as an integral over the element volume
where
where
3.2. Plane axi-antisymmetric finite element
An interesting application is represented by the study of axisymmetric structures subjected to axi-antisymmetric loadings (this term refers to the well known cases of symmetry and anti-symmetry). An example is a shaft of variable diameter under a torsion load applied at the ends (Timoshenko & Goodier, 1951). In this configuration, load is antisymmetric with respect to each plane passing across structure
In literature, to solve the problem of an axisymmetric body under torsion by using a plane FE approach it is suggested using a variational approach, in which the solution is a stress function φ which minimises the functional (Zienkiewicz & Cheung, 1965, 1967):
Instead, this work will show that it is possible to develop a simple FE theory applicable to an axisymmetric structure under torsion, by analogy with the theory of axisymmetric loading previously reviewed. In fact, in this configuration each element node has only one degree of freedom (the hoop displacement
where
where [
The expression of the element stiffness matrix and of equivalent nodal loads can thus be easily evaluated using Eq. (4), with the following elasticity matrix:
It is worth mentioning that this approach seems not to be implemented in commercial FE codes. Nevertheless, it can be implemented as a particular case of a harmonic finite element, discussed in next section.
3.3. Harmonic finite element
A third type of problem, of more practical interest, is when the structure is axially symmetric but the loading is not, so that the analysis is really three-dimensional. A great simplification can be obtained by using a semi-analytical approach, based on a harmonic FE model and Fourier series expansion of loads. As it will be shown, it can be demonstrated that, in linear analysis, a harmonic load produces a harmonic response in term of stress and displacements. The solution is then obtained by superimposing results of each harmonic, which are totally uncoupled (Cook et al., 1989; Zienkiewicz & Taylor, 2000).
To start with, the nodal loads applied to the structure can be expanded in Fourier series as:
where symbols
It is possible to demonstrate (Cook et al., 1989) that in a linear analysis, when loads are expanded as in Eq. (11), displacement components are described by Fourier series as well:
All three displacements are needed because the physical problem is three-dimensional. The motivation of the arbitrarily chosen negative sign in the
A Fourier series expansion similar to Eq. (12) can be equally used also for the nodal displacements of a finite element. Within a finite element, one can thus interpolate the amplitudes
where
The strain vector in cylindrical coordinate is given by (Cook et al., 1989):
where [∂] is a differential operator matrix, with dimension 6x3. Therefore, also strains are expanded in Fourier series and the contribution of
while for
Shape functions
Two stiffness matrices
The integrand matrix [
It should also be mentioned that, due to choice of negative sign in the second expression in Eq. (12), the stiffness matrix for double-barred terms is identical to that of single-barred terms, that is
where each matrix
4. Thermal analysis with a harmonic FE approach
An interesting issue here investigated is the use of a one-dimensional harmonic FE model for steady-state and transient thermal analysis of two-dimensional axisymmetric structures under non-axisymmetric thermal loadings. An example could be a rotating cylinder under imposed temperature and fluxes, which has been used as a simplified model for estimating the non-uniform transient temperature in a work roll of a hot rolling mill (Benasciutti, 2010a). Numerical FE modelling may be still very complex, even using a simplified two-dimensional modelling with a commercial FE code, as transient analysis needs large computational times and computer resources. Instead, a semi-analytical approach based on harmonic model would greatly reduce the computational burden. On the other hand, in commercial FE codes a one-dimensional harmonic thermal finite element is usually not implemented. Hence, this work will close the gap by developing the theory of a one-dimensional harmonic finite element for steady-state and transient thermal analysis of two-dimensional axisymmetric problems under non-axisymmetric thermal loadings. Other examples on thermal problems solved by a numerical approach can be found in (Awrejcewicz et al., 2007, 2009).
4.1. Steady-state thermal analysis
As discussed above, a harmonic model based on Fourier series expansion allows a three-dimensional physical problem to be reduced to a two-dimensional one. Similarly, a two-dimensional problem can be solved by a one-dimensional analysis. A significant reduction in total simulation time and also a saving of computational resources is thus achieved.
The study of two-dimensional problems by harmonic model needs a one-dimensional mesh along the radial direction. In a thermal analysis, the elements used are of "truss type" (one-dimensional) and the degree of freedom in each node is temperature. The theory of different types of harmonic finite element for thermal analysis have been formulated for a plane thermal analysis: two-node elements (with one or two Gauss points) having linear shape functions, three-node elements (with two or four Gauss points) having quadratic shape functions (Benasciutti et al., 2010b, 2011).
An example is shown in Fig. 3(a); two reference systems are used:
As previously discussed with reference to mechanical analysis, due to orthogonality property of trigonometric terms of Fourier series expansion the element stiffness matrix is a block diagonal matrix, see Eq. (17), where [
Explicit expressions for shape functions and stiffness matrix have been derived for each element type mentioned above (Benasciutti et al., 2011). As an example, the stiffness matrix of the two-node element with two Gauss points is here calculated. Similarly to Eq. (12), the temperature is first expanded in Fourier series as:
where
Within an element, amplitudes
where
(for more clarity, explicit dependence on variable
where
where the two coordinates
Theoretically, it is expected that increasing the node number would improve the numerical accuracy of results at the expense of higher computational cost. Element characterised by the highest rate of numerical convergence would allow the use of coarser meshes (with less number of elements), with a considerable decrease in number of equations to be solved and hence a significant reduction in required computational burden. Therefore, the choice of the most suitable element for a selected problem represents a crucial step in the analysis.
A test has been performed to compare the performance (accuracy and convergence rate) of different elements. A reference thermal problem is repeatedly solved by different elements, having various mesh densities. Considering that computational burden is approximately proportional to the number of equations to be solved and, in turn, to the number of nodes in the mesh, a comparable computational time would be roughly achieved by different meshes having the same number of nodes. The comparison then assumes that different elements have an equivalent mesh if the number of nodes is the same. For example, a mesh with 10 two-node elements (21 nodes) has approximately an equivalent computational burden to a mesh with 20 three-node elements (21 nodes).
The geometry and thermal load configuration used in test is shown in Fig. 4(a). A constant thermal flux is applied over the surface angular sector (π/8
The comparison is made with reference to the maximum temperature calculated by different element and mesh types. Figure 4(b) shows the trend as a function of the number of nodes in the mesh (a uniform element density is used); figures quoted in parenthesis indicate, as a rough guide, the computational time required by each element. The asymptotic value of temperature (for an infinitely fine mesh) has to be interpreted as the true (converged) value.
At a first sight, it appears clear how three-node elements give a much faster convergence rate, compared to two-node elements. The best performance (slowest running time) is given by three-node element with two Gauss points, which has been used in all subsequent simulations (see Fig. 9(a) ). Contrarily to what is expected, an increase of Gauss points (keeping the element nodes fixed) does not cause any increase in results accuracy (that is to say, in coarser meshes better results are obtained with lower Gauss points).
Another important parameter is the number of harmonics
Different types of imposed boundary conditions (temperature, flux and convection) have been studied in detail. Similarly, a thermal flux applied on boundary elements is converted into equivalent nodal loads. Instead, special attention has been deserved to some particular boundary conditions, as application of different temperatures on different boundary sectors or application of convective heat exchange, see Section 4.3.
4.2 Transient thermal analysis
The differential equation of equilibrium governing a transient thermal FE analysis is:
where
which depends on shape function matrix [
Explicit expressions for "mass matrix" have been calculated for two- and three-node elements (Benasciutti et al., 2011). For example, for a two-node element it is:
where
In a FE model, the "mass matrices" for each element in the mesh are assembled, to get a global "mass matrix" [
Mathematically, Eq. (23) represents a system of linear differential equations of second order, which in FE procedures are usually solved by numerical methods. In fact, in FE approach time domain is represented by a discrete sequence of time instants, in which solution is calculated. The time difference between adjacent time instants is the time step Δ
Numerical methods differ in the way the time derivative is approximated. The partial differential equation governing a transient thermal analysis can be solved by a finite difference method, implemented as an explicit or implicit numerical algorithm: explicit methods calculate the state of the system at a later time only considering the state of the system at current time, while implicit methods calculate the state of a system at a later time by solving an equation involving both the current state of the system and the later one. Implicit methods are usually slower (although more accurate) than explicit methods on single time step computation, as they require to solve a linear system at each time step. Conversely, explicit methods are usually faster, as they do not need to invert matrices.
In the present study, time integration has been performed by two different numerical methods. The first is
In FFD method, the time derivative in Eq. (23) is approximated as the discrete increment
Note that "mass matrix" [
An alternative and completely new original method has been developed for time integration. In analogy with the linear acceleration method (Bathe, 1996), the method here presented assumes a linear variation of the first derivative, thus it has been called
where
Unlike FFD algorithm, this method is implicit, as the system in Eq. (28) has to be solved at each time step. Furthermore, LSM is also "unconditionally stable" (i.e. solution converges independently of the choice of time step). Note that to further improve the computational speed, the system in Eq. (28) can be solved by LU factorization.
A test has been performed, to compare the relative performance and accuracy of FFD and LSM algorithm in semi-analytical approach. The geometry and thermal load configuration shown in Fig. 4(a) is analyzed under a physical transient of 900 s, starting from a uniform temperature of 0°C. The comparison is made both with a plane FE and a semi-analytical model (implemented with FFD method and LSM algorithm, with and without LU factorization), with
The same number of computation steps must be imposed in all FE modes (at least, in 2D and 1D-LSM), to allow an effective comparison of simulation time. A time step Δ
In harmonic model, both algorithms (FFD, LSM) have been shown to provide almost coincident results, which are also identical to those of 2D model. A time saving for LSM of about 4% compared to FFD and even 99% with respect to a plane FE model is observed. The use of LU factorization in LSM algorithm further reduced computation time of about one third. The above test then revealed that the fastest algorithm for transient thermal analysis is LSM method with LU factorization; this has been used in the illustrative example.
4.3. Boundary conditions
4.3.1. Prescribed thermal flux and temperature
In a thermal FE analysis, a prescribed temperature in a node is equivalent to an imposed nodal displacement in a structural analysis. In a harmonic model, the temperature value imposed on a node would be expanded by the Fourier series to the
This particular boundary condition, however, poses particular numerical problems, due to one-dimensional nature of harmonic model. Three methods have been proposed (Benasciutti et al., 2011): two of them (fictitious thermal flux, imposed temperature) are based on Fourier series expansion of the imposed temperature, the third is based on Lagrange multipliers.
In the first method, the applied temperature is interpreted as a consequence of a fictitious thermal flux (unknown) on the same surface portion, which has to be determined. The prescribed temperature is expanded in Fourier series, similarly to thermal flux; the coefficients for thermal flux expansion are then calculated from those defining the series expansion of applied temperature. In the second method, the prescribed temperature is treated similarly to an imposed displacement (constraint) in a mechanical analysis, in which the equation corresponding to the prescribed degree of freedom is cancelled out in the system of equilibrium equations. In the harmonic model, this procedure has to be applied to every term in Fourier series expansion. In the third method with Lagrange multipliers, the "stiffness matrix" is calculated by a variational approach, which minimises the total potential energy Π. In stiffness matrix [
where
Illustrative examples showed that the performance (accuracy and computation speed) of all three methods is quite comparable, although the third one gives more flexibility in representing various boundary conditions. In particular, an imposed temperature and flux over the same boundary portion is not possible with the first two methods. The decisive advantage of the method using Lagrange multipliers is its possibility to represent any variation with angle
4.3.2. Convection
A convective flux per unit area is defined as
Three different methods with various approximation levels have been developed: one based on constant convective coefficient and single average surface temperature, one based on a step-wise constant convective flux and average surface temperature, the third based on trigonometric formulae (which has been implemented by two different algorithms: temperature calculated either at previous or at current time step) (Benasciutti et al., 2011).
In the first method, the constant convective coefficient α over an angular sector (
where
The third method tries to capture the real surface temperature variation within sector (
The product of the two series is further simplified through the well-known prosthaphaeresis trigonometric formulae, arriving at the following expression (for node
where [Αi] is a matrix including amplitudes of Fourier series of α (it can be then calculated only once) and
where {
5. Numerical examples
The numerical examples here discussed address both a mechanical and thermal analysis, with the aim of testing the performance and accuracy of semi-analytical approach described above. Two different applications are discussed: a shaft under torsion, bending and axial mechanical loadings and a rotating cylinder under surface thermal loadings. The study confirms how semi-analytical approach is able to provide accurate results, with in addition a strong reduction in computation time, compared with classical two- or three-dimensional FE simulation.
5.1. Mechanical analysis
A first case study here presented refers to a cylindrical component with a shoulder fillet, loaded under torsion, bending and axial load applied at the ends, see Fig. 6(a). The geometry shown is a type of stress concentration frequently encountered in shafts, axle spindles and rotors machine design. The shoulder fillet is designed to have
This geometry can be solved by a harmonic model; the axial, bending and torsion loads can be represented by surface loads expanded in Fourier series as:
where
The axisymmetric structure under the three different loading configurations can be easily solved by using three plane models, see Fig. 6(b)-(d): an axisymmetric model for the axial load and an axi-antisymmetric model for torsion; bending can be treated by means of a harmonic model, with only one term in Fourier series (in fact, bending stresses follow a sinusoidal variation law with respect to symmetry axis).
The figure shows the plane mesh used, as well as the applied loads and constraints. The adequacy of the mesh has been tested with a convergence analysis based on the comparison of stress concentration factors calculated by increasing levels of mesh refinement. The numerical solution given by a three-dimensional FE model with hexahedral elements is also used for comparison purposes. The stress concentration factors numerically calculated by plane models, shown in Fig. 6, have been compared with theoretical values from manuals (Peterson, 1974), as well as with numerical values from three-dimensional FE analysis. The overall comparison shows a general good agreement. Compared to three-dimensional model, the semi-analytical approach greatly reduces the overall computation time.
The semi-analytical approach here described can thus be extended to the analysis of a shaft loaded axially, in torsion and in bending, showing that the stress distribution can be accurately described (at least far from the points where loadings are applied) by using only three terms of the Fourier series, i.e. the constant and the first harmonic term.
A second case study here analysed is thus a simply supported shaft, loaded by non-planar distributed loads, see Fig. 7. A three-dimensional model and a plane harmonic model, both having similar mesh distributions, have been considered.
Table 1 shows the maximum tensile stress
localised re-meshing | no localised re-meshing | |
30 | 1.00 ( | 0.92 ( |
7 | 1.00 ( | 0.92 ( |
3 | 0.98 ( | 0.91 ( |
3D model | 1.00 ( | 0.84 ( |
Results are reported respectively for the original mesh and after a local refinement close to stress concentrations. It is possible to observe that the plane model with only 2 harmonics and non-refined mesh gives better results with respect to those achievable with the three-dimensional model with the same mesh distribution. In the case of a plane model with localised re-meshing, a relevant reduction in the error can be observed, without a significant increment of computational time, even if the number of degrees of freedom increases of about 10 %. It can be noticed that, with only 3 terms of Fourier series, an error lower than 2% can be obtained. Convergence can thus be achieved quite easily with the harmonic model. With respect to the harmonic case, the three-dimensional model allows to achieve a similar accuracy only at the expense of unfeasible computational times.
5.2. Thermal analysis
The second example refers to transient analysis of a rotating cylinder under thermal loadings (see Fig. 8), which has been used as a simplified model for a work roll of a hot rolling mill (Benasciutti et al., 2010a). In the simplified approach, only work roll (without the strip) was modelled in a two-dimensional analysis; thermal loadings are applied on the surface, to simulate strip heating and water jet cooling.
The thermal load configuration used in simulations is shown in Fig. 8(a): an infinitely-long cylinder, rotating at constant angular speed and subjected to constant input heat flux and convective cooling. The numerical analysis simulates a thermal transient of 3600 seconds, which corresponds to about 1690 roll revolutions. At the initial simulation time, work roll is assumed at a constant uniform temperature
The work roll configuration in Fig. 8 is an example of axisymmetric structure under non-axially symmetric thermal loadings, which can be solved by the semi-analytical approach previously described. A two-dimensional FE model, implemented by the commercial code ANSYS®, is also used for comparison purposes, to test the accuracy of harmonic FE model.
Figure 8(b) shows the plane FE model of work roll. A mesh refinement is imposed near the surface, along the tangential and radial directions, to capture the thermal gradient here expected. Small elements are located for a depth of 10% of work roll radius, with even smaller elements placed immediately underneath the surface, for a depth of 2% of radius.
Since in FE analysis rigid body motion is not allowed, work roll rotation has been simulated by considering the roll at rest and by applying rotating thermal loadings. As an order of magnitude, the simulation required about 3 days of simulation. It is worth mentioning how results of thermal analysis have been validated, see (Benasciutti et al., 2010b), by an analytical solution for the stationary temperature distribution (Patula, 1981). Other details on FE model and simulation parameters can be found in (Benasciutti 2010a).
Cylinder radius | |
ω=2,953 rad/s | Cylinder angular speed |
Heating sector | |
β = 45° | Angular gap between heating and cooling |
Ψ = 90° | Cooling sector |
Input thermal flux | |
α=10100 W/m2K | Convection coefficient |
Bulk temperature of cooling medium | |
Initial temperature of work roll |
In semi-analytical approach, the 2D geometry of work roll in Fig. 8(a) is represented by a 1D model. A three-node element with two Gauss points is used, see Fig. 9(a); the mesh, shown in Fig. 9 (b), has a total amount of 28 elements and 57 nodes. Note that, for an effective comparison, both 1D and 2D models have an identical element distribution in radial direction, compare Fig. 8(b) with Fig. 9 (b). In addition, the same time step in transient analysis has been chosen for both models, to allow a comparison of running times.
In harmonic model, the choice of the correct number of harmonic
For simulation of convective cooling, the choice is between two algorithms based on trigonometric functions. For
Conversely to plane FE approach, which needs that work roll must be fixed and thermal loads rotating, in harmonic model two options are available: work roll fixed with rotating thermal loading, or vice versa. A comparative study (Benasciutti et al., 2011) showed that both methods give similar running times, with the first one (work roll fixed) slightly faster (and also more simple to implement) than the second. Therefore, method with rotating thermal loading has been preferred.
Results for thermal simulations are presented from Fig. 10 to Fig. 13. For example, Fig. 10 shows the temperature field in work roll after 1800 seconds calculated by plane FE model. The temperature map at other time instants (not included here) would show a progressive heating of the entire work roll, with the largest temperature gradients localised very close to the surface (which justifies the use there of very small elements in the mesh). The temperature field calculated by semi-analytical FE model is not provided, as it is very similar to that given by plane FE model.
Figure 11, instead, compares the temperature time history within a 60-seconds time interval, for a point located on work roll surface. Each peak temperature occurs when the monitored point enters the heating zone, thus the series of equally-spaced peaks identifies the sequence of work roll rotations. The progressive increase of peak temperature in consecutive rotations confirms the transient nature of the thermal phenomenon here investigated. A very similar trend is observed for both FE models; the small differences may be attributed to the different rate of results saving on computer hard disk.
The result shown in Fig. 12 refers instead to the radial temperature distribution at next time instants. Only temperature distribution for plane FE model is shown, as that calculated by harmonic model would be practically coincident. The continuous temperature increase with time, especially inside the work roll, is indicated by the different curves. Both diagrams also
confirm that work roll remains at a rough uniform temperature, except for a very small portion very close to the surface (wide about 1% of roll radius), where a steep temperature gradient is observed (note that diagrams only plot radial coordinates close to roll surface). The extremely localized nature of temperature variation within this narrow region is usually called "thermal boundary layer" in technical literature.
Finally, Fig. 13 shows the temperature change for points at different radial depths, along a fixed angular position. A close agreement between semi-analytical and plane FE models is confirmed. The figure further emphasizes that the thermal gradient is confined in a small region close to the boundary. In fact, on work roll surface the temperature ranges of about 200°C, while at 1 mm depth from surface the variation is only 100 °C, and even negligible at 6 mm below surface. In addition, the detail in the same figure highlights a sort of "thermal inversion" phenomenon induced by forced convection cooling, in which the work roll material on the surface is at a lower temperature than material inside. The lateral expansion of surface elements, prevented by the surrounding elements at lower temperature, is the basic mechanism which explains the development of thermal stresses in work roll.
6. Conclusions
This work has investigated theory and application of simplified FE approaches for the analysis of axisymmetric structures loaded by non-axisymmetric loadings. The aim was to develop alternative FE methods, which allow obtaining the solution of complex three-dimensional problems through a combination of several simpler and faster one- and two-dimensional analyses, which usually require reduced computational efforts. The study has focused on both mechanical and thermal problems, in which the structure is axisymmetric, but the load is not.
In the mechanical context, the classical axisymmetric stress analysis has been first reviewed from literature, as it constitutes a reference example for the subsequent discussion. Using a similar approach, an original finite element is next developed, for solving axisymmetric structures under axi-antisymmetric applied loads, as for example a shaft under a torsion load. General equations for displacement field, strain and stress have been derived.
Another and more general class of problems is that of axisymmetric structures loaded non-axisymmetrically, which are solved by a semi-analytical finite element approach based on Fourier series expansion of applied loads and displacement field. In a linear analysis, the harmonics in Fourier series are totally uncoupled due to orthogonality property of trigonometric functions. Therefore, a complex three-dimensional problem is solved by superposition of solutions of several two-dimensional analyses, or similarly a two-dimensional problem is replaced by several one-dimensional analyses. A considerable reduction in computational times and computer resources is then achieved. The theoretical framework for harmonic finite element analysis has been first presented, as the fundamental equation for stiffness matrix and the stress/strain relations. The performance of the semi-analytical approach is discussed with a numerical example: a shaft with a shoulder filled under respectively an axial, bending and torsion loading. The stress concentration factors calculated by harmonic model are compared with values provided by manuals. The study has confirmed that harmonic model is capable to predict with great accuracy the stress concentration factors of a three-dimensional geometry under three different loading conditions, by using simple two-dimensional models. A second case study showed that the case of a shaft loaded by non-planar loads can be solved quite easily by using a harmonic model with only three terms of the Fourier series. This approach gives a relevant advantage in terms of computational time compared to three-dimensional modelling.
For what concerns thermal problems, it has been developed a harmonic finite element approach for the steady-state and transient thermal analysis of two-dimensional axisymmetric structures under non-axisymmetric thermal loadings. The theory of different types of one-dimensional finite element has been derived: two-node elements (with one or two Gauss points) or three-node elements (with two or four Gauss points). The relative performance of all elements has been compared in terms of accuracy and computation speed; the best performance is provided by three-node element with two Gauss points. Also the choice of the number of harmonic is investigated with a benchmark test.
For transient thermal analysis, explicit expression of "mass matrix" for the two-node element has been presented. Furthermore, two algorithms for numerical time integration has been introduced: the
The presented results have confirmed the great reliability offered by semi-analytical approach, in providing accurate results with at the same time a significant reduction in computational times, compared to classical FE analyses.
References
- 1.
Awrejcewicz J. Krysko V. A. Krysko A. V. 2007 Thermo-dynamics of plates and shells Springer-Verlag,3-54034-261-3 Germany). - 2.
Awrejcewicz J. Pyryev Yu. 2009 Nonsmooth dynamics of contacting thermoelastic bodies Springer Science,978-0-38709-652-0 New York (USA). - 3.
Bathe K. J. 1996 Finite element procedures in engineering analysis Prentice-Hall,0-13301-458-4 Jersey (USA). - 4.
Benasciutti D. Brusa E. Bazzaro G. 2010 Finite element prediction of thermal stresses in work roll of hot rolling mills. 2 1 707 716 - 5.
Benasciutti D. De Bona F. Munteanu M. Gh. 2010b Harmonic model for numerical simulation of thermal stresses in work roll of hot rolling mill. , Paris, May 2010. - 6.
Benasciutti D. De Bona F. Munteanu M. Gh. 2011 Work roll in hot strip rolling: a semi-analytical approach for estimating temperatures and thermal stresses. To be included in , Mali Losinj, Croatia, June, 2011. - 7.
Bressan F. De Bona F. Munteanu M. Gh. 2009 Semi-analytical Finite element for Shaft Design, , Dubrovnik, Croatia, September-October 2009. - 8.
Cook R. D. Malkus D. S. Plesha M. E. 1989 (3rd ed.), Wiley & Sons,0-47184-788-7 - 9.
Genta G. Tonoli A. 1996 An harmonic finite element for the analysis of flexural, torsional and axial rotordynamic behaviour of discs. ,196 1 19 43 - 10.
Kim J. R. Kim S. J. Kim W. D. 1994 Parallel computing using semi-analytical finite element method. ,32 5 1066 1071 - 11.
Lai J. Y. Booker J. R. 1991 Application of discrete Fourier series to the finite element stress analysis of axi-symmetric solids ,31 4 619 647 - 12.
Patula E. J. 1981 Steady-state temperature distribution in a rotating roll subject to surface heat fluxes and convective cooling ,103 36 41 - 13.
Pedersen P. Laursen C. L. 1982 Design of minimum stress concentration by finite element and linear programming. J Struct. Mech.,10 4 375 391 - 14.
Peterson R. E. 1974 . Wiley & Sons,0-47168-329-9 York. - 15.
Taiebait H. A. Carter J. P. 2001 A semi-analytical finite element method for three-dimensional consolidation analysis. ,28 1 55 78 - 16.
Thomas T. J. Nair S. Garg V. K. 1983 Elasto-plastic stress analysis and fatigue life prediction of a freight car wheel under mechanical and cyclic thermal loads ,17 3 313 320 - 17.
Timoshenko S. P. Goodier J. N. 1951 (2nd ed.), McGraw-Hill, USA. - 18.
Wilson E. L. 1965 Structural analysis of axisymmetric solids. ,3 12 2269 2274 - 19.
Zienkiewicz O. C. Cheung Y. K. 1965 Finite element in the solution of field problems. ,24 507 510 - 20.
Zienkiewicz O. C. Cheung Y. K. 1967 Stresses in shaft. ,24 696 697 - 21.
Zienkiewicz O. C. Taylor R. L. 2000 (5th ed.), Butterworth-Heinemann,0-75065-049-4