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 *θ* of a cylindrical reference frame (*r*,*θ*,*z*), where *z* is the component axis of symmetry and *r* is the radial distance from *z*-axis, see Fig. 1.

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 *z*-axis, the analysis is axisymmetric and mathematically two-dimensional, that is results are independent of *θ* and they are only function of *r*, *z* coordinates (Bathe, 1996; Cook et al., 1989; Zienkiewicz & Taylor, 2000). Examples are disks rotating at uniform speed under centrifugal forces or cylindrical vessels under internal pressure.

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 *z*-axis and also independent of *θ* angle. An example is a cylindrical body under a torsion loading. As it will be shown, the analysis becomes really one-dimensional and results are only function of *r* and *z*.

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 *u* and *w* in any plane section of the body along its axis of symmetry completely define the state of strain and, accordingly, the state of stress (Zienkiewicz & Taylor, 2000). Thus, the circumferential (hoop) displacement *v*, the tangential stress components *τ*_{rθ} and *τ*_{θz} and their corresponding shear strains γ_{rθ} and γ_{θz} must be zero, see Fig. 2. The analysis then reduces to a plane FE model, characterized by only radial *u*(*r*,*z*) and axial *w*(*r*,*z*) displacements, where *r* and *z* denote the radial and axial coordinates of a point within the structure. For a *M*-nodes finite element, the vector of displacement field in the cylindrical reference system (*r*,*θ*,*z*) is:

where *u*_{i}, *w*_{i} are the nodal displacements and *N*_{i}(*r*,*z*) is the shape function for node *i*.

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 [*B*_{i,pl_st}] corresponds to that of plane case and its dimension is 3x2. For example, as it is well known for 3-nodes triangular element it is:

where *A*=(*r*_{1}*z*_{3}+*r*_{2}*z*_{1}+*r*_{3}*z*_{2}–*r*_{1}*z*_{2}–*r*_{2}*z*_{3}–*r*_{3}*z*_{1})/2 is the element area. As matrix [*B*] depends on *r* and *z* coordinates, the strains are no longer constant within an element, as in plane stress or plane strain (the strain variation is due to the *ε*_{θ} term).

The element stiffness matrix is calculated as an integral over the element volume *V*_{el}, which for axial symmetry coincides with the whole ring of material:

where *A*_{el} is the cross-sectional area of the element on a plane section. The elasticity or Hookean matrix [*D*] in Eq. (4), which links the vectors of strains and stresses, in the hypothesis of isotropic material has the following form:

where *E* is the modulus of elasticity (Young's modulus) and ν is the Poisson's ratio.

### 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 *z*-axis (this plane then behaves as a plane of axis-antisymmetry) and it is also independent of *θ* angle.

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 *v*), while radial and axial displacements *u* and *w* (warping), as well as normal stresses *σ*_{r}, *σ*_{θ}, *σ*_{z}, shear stress *τ*_{rz} and their related strain components vanish. By symmetry, the hoop displacement does not depend on angle *θ* and only two non-null strains *γ*_{rθ} and *γ*_{θz} are present*.* By analogy with Eq. (1), the displacement of a point within a *M*-nodes element is:

where *v*_{i} is the nodal displacement and *N*_{i}(*r*,*z*) is the shape function of node *i*. The related strain vector is:

where [*B*] is the strain-displacement matrix, with dimension 1x*M*. As an example, for a 3-nodes triangular finite element, matrix [*B*] has the following expression:

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 *R*, *T* and *Z* indicate respectively the radial, hoop and axial load components. A similar series expansion holds also for body forces, boundary conditions, initial strains, etc. In Eq. (11), all barred quantities are amplitudes, which are functions of *r*, *z* but not of *θ*. Single-barred amplitudes represent symmetric load components (loads which have *θ*=0 as a plane of symmetry), while double-barred amplitudes represent antisymmetric load terms. The sine expansion in *T* load is necessary to assure symmetry, as the direction of *T* has to change for *θ*>π. The constant terms *R*_{0} and *Z*_{0} permit axisymmetric load condition to be described, while the term *T*_{0} refers to the axi-antisymmetric load. These three terms are grouped into a single vector representing the constant term of the Fourier series.

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 *v* series is that it greatly simplifies the computation of the element stiffness matrix, as it will be explained later on. As for the loads, the single- and double-barred terms refer to symmetric and antisymmetric components.

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*in* specifies that amplitude refers to node *i* and harmonic *n*. Therefore, the vector of displacement field within the element can be described in the following form:

where *M* arrays of 3x3 submatrices (*M* is the number of element nodes); note that such matrices both depend on *n* because of the (cos*nθ*) and (sin*nθ*) terms.

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 *n*-th harmonic thus is*n*, the strain-displacement matrices *M* arrays of the 6x3 submatrices *i*-th node, one submatrix is:

while for *nθ*) and (cos*nθ*) are interchanged and, in addition, there is an algebraic sign change in the third and fourth row, that is to say that *nθ*) with (cos*nθ*) and (cos*nθ*) with (sin*nθ*).

Shape functions *N*_{i} depend on *r* and *z*. Therefore, as shown in Eq. (15)*r*, *z*, *θ* and obviously of *n*. A unique strain-displacement matrix can be defined by assembling matrices *H* harmonics are retained in the Fourier series, matrices *H* arrays of the 6x3*M* submatrices

Two stiffness matrices

The integrand matrix [*B*]^{T}[*D*][*B*] is a full matrix of size (3*MH*)x(3*MH*); it is composed of an *H*x*H* array of 3*M*x3*M* submatrices [*k*_{n}]. The off-diagonal submatrices contain in every term the products (sin*mθ*)(sin*nθ*) or (cos*mθ*)(cos*nθ*) with *m*≠*n*, which give zero when integrated from –π to +π, due to the so-called orthogonality property of trigonometric functions. The remaining *H* on-diagonal submatrices, with dimension 3*M*x3*M*, contain (sin^{2}*nθ*) and (cos^{2}*nθ*) in every term, which integrated from –π to +π give a common factor π (or 2π for *n*=0). Integration on *r* and *z* variables in Eq. (16) is done, as if the problem were axially symmetric.

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 *M*x3*M* and subscript 0,1, …,*H* specifies the number of the Fourier harmonic. Matrix [*k*_{0}] has similar dimension; it contains both symmetric and antisymmetric terms, but obviously no coupling terms are present. The above theory then shows that in linearity, due to the topology of stiffness matrix, the problem is uncoupled and *H* separate problems are solved (Bressan et al., 2009). For example, if the simple case of a cylinder is considered, in a three-dimensional analysis it would be necessary to solve about (3*p*^{3/2})/(*a*^{1/2}) equations (parameter *a* is the cylinder aspect ratio, defined as the ratio between the length and the radius, parameter *p* is the plane model node number), whereas in the plane case (2*H*+1) systems constituted by 3*p* equations should be solved. As *H*<<*p*, computational time would be strongly reduced.

## 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: *r* is the abscissa in global reference system (*r*=0 is the centre of the axisymmetric geometry), *x* is the coordinate within the element. The mesh is one-dimensional and consists of adjacent elements located along the radius of solid, see Fig. 3(b).

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 [*k*_{n}] is the elementary stiffness matrix for *n*-th term of Fourier series and *H* is the number of harmonics.

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 *r*, but not on *θ*. In practice, only a finite number of harmonics *H* is used in the summation in Eq. (18).

Within an element, amplitudes *nodal* amplitudes

where *N*_{1}(*r*) and *N*_{2}(*r*) are the shape functions, which are linear in a two-node element: *N*_{1}(*r*)=(*r*_{2}–*r*)/*L* and *N*_{2}(*r*)=(*r*–*r*_{1})/*L*. After some simple matrix algebra, the vector of "strains" (derivatives of temperature) in polar coordinates can be written in matrix notation as:

(for more clarity, explicit dependence on variable *r* is omitted). Element stiffness matrix for harmonic *n* is calculated as in Eq. (4) or (16). Matrix product inside the integral contains terms as (sin^{2}*nθ*) and (cos^{2}*nθ*), which integrated from –π to +π give a factor π (or 2π for *n*=0). The following expression is then obtained (Benasciutti et al., 2011):

where *h* is element thickness in *z* (axial) direction, [*D*]=diag(*λ*, *λ*) is a diagonal matrix with material thermal conductivity *λ*. The integral in Eq. (21) can be solved numerically (Gauss quadrature), obtaining a closed form solution. For example, for the two-node element it is:

where the two coordinates *x*_{G} and *r*_{G} are used to specify the position of Gauss point in the element and global reference system, respectively. As it can be seen, element stiffness matrix depends on both harmonic *n* and also on radial position of Gauss point, *r*_{G}. This implies that elements located at different radial positions have different stiffness matrix. Analogous expressions can be easily obtained for other elements, similarly to what done above. The increase of Gauss points would improve the accuracy of numerical integration in Eq. (21).

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*H*=12 harmonics is used in Fourier series expansion.

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 *H*, which should be chosen as the best balance between accuracy and simulation time. Higher *H* values are expected to give higher precision in Fourier expansion of applied loads (especially for step loads), as well as better accuracy and higher simulation time. A test is performed to assess the effect of the number of harmonic *H* on accuracy of results and computational speed. Results from semi-analytical model with different discrete values of *H* (i.e. 4, 8, 12, 18, etc.) have been compared with those by a two-dimensional FE model (commercial code ANSYS^{®} has been used). A relative error between the maximum temperature calculated by harmonic model and that of plane FE model is also calculated. The comparison is shown in Fig. 4(c); with *H*=12 the absolute relative error is about 1%, as can be seen, while for higher harmonics the error rapidly decreases below 0.2%. On the other hand, it has been noted that simulation time in harmonic model increases roughly linearly with *H* (because equilibrium equations are not coupled). Also the approximation error (Gibbs' phenomenon) in Fourier series expansion of step-wise periodic functions should be considered carefully in selecting the optimal value of *H*. For the example examined, the best choice would be *H*=24, which, however, does not actually represent a result of general validity. In fact, each problem needs to be carefully analyzed to find the best *H* value.

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 *M*] and [*K*] are the "mass" and "stiffness" matrices (assembled from the corresponding element matrices), {*F*} is the vector of externally applied loads. For a one-dimensional element, the "mass matrix" is defined as:

which depends on shape function matrix [*N*(x)]=[*N*_{1}(x) *N*_{2}(x)], as well as on material volumic mass *ρ* and specific heat *c*. Unlike matrix [*B*_{n}], the "mass matrix" does not explicitly depend on the harmonic *n*, hence it can be calculated only once for all terms of Fourier series expansion.

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 *r*_{1} is radial position of node "1" in global reference system. As shown, [*m*_{el}] depends on the position of element in absolute reference system, which means that two elements of equal length, but different location in the mesh, would have a different "mass matrix".

In a FE model, the "mass matrices" for each element in the mesh are assembled, to get a global "mass matrix" [*M*], similarly to assemblage of global "stiffness matrix" [*K*]. Similarly to [*K*], also [*M*] is a block diagonal matrix, with non-zero terms close to the main diagonal.

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 Δ*t*.

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 *Forward Finite Difference* (FFD) method, which is an explicit method proven to be very quick and effective. However, it is said to be "conditionally stable", as it requires the time step Δ*t* be smaller than a critical value to get a stable solution. This can represent a serious disadvantage, as the need of a stable solution may require a relatively small time step, which can result in a very large total simulation time.

In FFD method, the time derivative in Eq. (23) is approximated as the discrete increment*u*}_{i} and {*u*}_{i+1} are the vectors of nodal temperatures calculated at consecutive time steps *i* and *i*+1. Substituting in Eq. (23) and rearranging, it follows the fundamental equation for FFD algorithm:

Note that "mass matrix" [*M*] is time-independent, thus it can be inverted only once, with a considerable time saving.

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 *Linear Speed Method* (LSM); in symbols:

where *i* and *i*+1, while

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 *H*=24. Additional material parameters used in simulation are volumic mass (*ρ*=7800 kg/m^{3}) and specific heat (*c*=460.5 J/kgK). Plane model adopts eight-node thermal elements (which have 3 nodes on each side), while harmonic model uses 12 elements with three nodes and two Gauss points (therefore, both models have the same number of nodes along the radial direction).

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 Δ*t*=1 s is chosen, so that the interval of 15 minutes (900 s) is divided into 900 time steps. Unfortunately, such a value is shown to be inadequate for 1D-FFD model; it is greater than the critical time step, thus it causes the numerical solution to not converge (numerical instability). A lower time step Δ*t* must then be imposed; for example, an acceptable value is Δ*t*=0.0335 s, although it gives a much greater number of time steps (26866) for the same physical transient of 15 minutes here analyzed.

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 *whole* circumference (nodal circle) passing through the node. Thus, it appears that in harmonic model only a constant temperature could be applied on the boundary, while in practical application it is often necessary to prescribe different temperature values along the boundary.

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 [*k*], formed by submatrices [*k*_{ni}], the generic element would be *k*^{ij}=∂^{2}Π/∂*u*^{i}∂*u*^{j}, where indexes *i*, *j* range from 1 to *nF*∙*M*, with *M* total number of nodes and *nF* number of Fourier terms. Without going into details, it suffices to say that the method seeks the minimum of a total potential energy

where *β* is the Lagrange multiplier. As a result, the stiffness matrix calculated by this third method is no longer uncoupled, as additional new rows and columns are inserted to describe the imposed temperature.

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 *θ* for the prescribed temperature.

#### 4.3.2. Convection

A convective flux per unit area is defined as *q*_{α}=α(*u–u*_{∞}), where α is convection coefficient, *u* the surface temperature and *u*_{∞} a reference (bulk) temperature of surrounding medium. A convective flux then depends on the surface temperature and may also change over time. Then, algorithms for transient analysis can be used to solve convective boundary condition.

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 (*θ*_{1}_{2}) is expanded in Fourier series, similarly to an imposed thermal flux. The convective heat flux is:

where *u*_{av} is the average surface temperature (calculated at previous time iteration), *R* is the radius of the axisymmetric plane body and *θ*_{1}*θ*_{2}); obviously, accuracy decreases as surface temperature shows large variations over interval (*θ*_{1}*θ*_{2}). An improvement can be obtained by further dividing (*θ*_{1}*θ*_{2}) into sub-sectors; a step-wise constant convective coefficient and average temperature are then calculated in every sub-sector.

The third method tries to capture the real surface temperature variation within sector (*θ*_{1}*θ*_{2}). Two separate Fourier series are then used, one describing surface temperature relative to bulk temperature *u*_{∞}, one modelling convective coefficient:

The product of the two series is further simplified through the well-known prosthaphaeresis trigonometric formulae, arriving at the following expression (for node *i*):

where [Α^{i}] is a matrix including amplitudes of Fourier series of α (it can be then calculated only once) and *u*. Vector *A*^{i}] are assembled to get a global sparse matrix [*Α*], with same dimension as [*K*] and [*M*]. Equation (23) for transient thermal analysis can be updated with convective flux:

where {*u*} can be calculated either at current or next time step. Without going into the details, it suffices to say that a numerical test was performed to compare the relative performance (accuracy and speed) of methods. The method with constant convective coefficient was shown to be not very accurate, as it tends to over-estimate the actual flux where surface temperature is lower. A slight improvement (at the expense of higher computational time) is achieved with a step-wise constant flux, although an increase in the number of surface divisions may cause numerical instability (a smaller Δ*t* must then be used). However, the best accuracy is reached by the method with trigonometric functions, which also has the lowest computational time (which is slightly lower to that of first method). In addition, the test example used did not reveal any difference between the algorithms which computes temperature at current or next time step.

## 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 *D*/*d*=1.5, *r*/*d*=0.1, so that stress concentration factors given by literature (Peterson, 1974) are respectively *K*_{t,N}=1.88 (axial load), *K*_{t,Mf}=1.72 (bending), *K*_{t,Mt}=1.37 (torsion).

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 _{0,axial} and σ_{0,bend} represent the maximum stress for axial and bending loading, while τ_{0} is the maximum shear stress for torsion load.

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 *σ*_{z} close to a fillet radius, normalised to the asymptotic stress value *σ*_{z∞} calculated by a three-dimensional FE model with a very refined mesh. Figures in parenthesis indicate a rough estimate of computational time required by each analysis.

Terms of Fourier series | stress ratio, σ_{z}/σ_{z∞} (computational time) | |

localised re-meshing | no localised re-meshing | |

30 | 1.00 (90 s) | 0.92 (90 s) |

7 | 1.00 (21 s) | 0.92 (21 s) |

3 | 0.98 (9 s) | 0.91 (9 s) |

3D model | 1.00 (days) | 0.84 (12600 s) |

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 *T*_{roll}=20°C. Simulation parameters, assumed for simplicity as temperature invariant, are summarized in Table 2.

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).

Parameter | Value |

R=300 mm | Cylinder radius |

ω=2,953 rad/s | Cylinder angular speed |

ϕ = 10° | Heating sector |

β = 45° | Angular gap between heating and cooling |

Ψ = 90° | Cooling sector |

q_{0}=13,7×10^{6} W/m^{2} | Input thermal flux |

α=10100 W/m^{2}K | Convection coefficient |

T_{0}=20°C | Bulk temperature of cooling medium |

T_{roll}=20°C | 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 *H* is critical, especially for the load configuration in Fig. 8 (a), where a stepped flux is applied over a very small angle (=10°). In fact, in Fourier series expansion the approximation error tends to increase as step width decreases, because of Gibbs' phenomenon at jump discontinuity; the error could be minimized by an appropriate high number of harmonics. Several benchmark tests were performed (Benasciutti et al., 2011) to identify the optimal number of harmonics for the configuration here analyzed. A comparison of the maximum transient temperature, calculated by plane model and various harmonic models with different number of harmonics revealed that *H*=100 would be an optimal compromise between accuracy and computing time.

For simulation of convective cooling, the choice is between two algorithms based on trigonometric functions. For *H*=100 harmonics, the test revealed that algorithm based on temperature computed at previous time step is faster (of about 20%) than algorithm which computes temperature at current time step, although both give comparable accuracy levels.

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 *Fast Forward Difference* (FFD) and the completely original *Linear Speed Method* (LSM), also implemented by LU factorization. A benchmark test revealed that the fastest algorithm for transient analysis is LSM with LU factorization. An illustrative example is finally discussed, which refers to a rotating cylinder under thermal loadings, used as a simplified model for work roll in hot rolling mill. The transient temperature distribution is estimated by a semi-analytical approach and a plane FE model. The harmonic model has been shown to give results in very close agreement with those of plane FE model, with however a significant reduction of total simulation time.

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.