A New Pulley Stress Analysis Method Based on Modified Transfer Matrix
Xiangjun Qiu and Vinit Sethi, USA
Summary
A new pulley stress analysis method is presented. It shall be referred to as the Modified Transfer Matrix (MTM) method. This method is based on a reformulation of transfer matrices for the pulleys cylindindrical shell, enddisk plate with nonuniform thickness and shaft by using finite element concepts. It combines the strength of both classical stress analysis methods and finite element methods. It proves to be an efficient and effective approach in determining the stresses in a pulley. A pulley stress analysis software program named PSTRESS 3.0 has been developed based on this new method. At the end of the paper, a numerical example of the pulley stress analysis, using PSTRESS 3.0, is given. The result is satisfactorily compared with that obtained in a finite element model (ANSYS) solution with a very fine mesh.
Nomenclature
x Cartesian coordinate in horizontal direction y Cartesian coordinate in vertical direction z Cartesian coordinate in pulley axial direction r r cylindrical coordinate in pulley radial direction cylindrical coordinate in pulley circumferential direction P shear force acting on shaft cross section M bending moment acting on shaft cross section q distributed load acting on shaft axis W_{s} transverse displacement of shaft neutral axis θ_{s} rotational angle of shaft cross section EI shaft bending stiffness a shaft cross section area G shear modulus of shaft material β shaft state variable vector A matrix containing coefficients of governing ODEs for
shaft bending deformationB vector representing the nonhomogeneous term of the governing ODEs for shaft bending deformation T transfer matrix for the governing ODEs for shaft
bending deformationI identity matrix K_{BM} TMB beam element stiffness matrix U_{BM} TMB beam element displacement vector F^{BM}_{ext} TMB beam element external force vector F^{BM}_{int} TMB beam element internal force vector Q_{r} transverse shear force acting on the disk cross section
perpendicular to radial directionQ_{} transverse shear force acting on the disk cross section
perpendicular to circumferential directionM_{r} bending moment acting on the disk cross section
perpendicular to disk radial directionM_{} bending moment acting on the disk cross section
perpendicular to disk circumferential directionM_{r} twisting moment acting on the disk cross sections perpendicular
to disk radial direction and disk circumferential directiong external transverse force acting on the neutral surface of the disk u pulley disk and shell displacements in pulley axial direction v pulley disk and shell displacements in pulley circumferential direction w pulley disk and shell displacements in pulley radial direction E YOUNG'S modulus of disk or cylindrical shell of pulley POISSON'S ratio of disk or cylindrical shell of pulley t thickness of disk or cylindrical shell of pulley c coefficient to describe the geometry of disk variable thickness p exponential number to describe the geometry of disk variable
thicknessD cross section bending stiffness of disk or cylindrical shell m FOURIER component number u_{m} FOURIER component of u v_{m} FOURIER component of v w_{m} FOURIER component of w g_{m} FOURIER component of g Q_{rm} FOURIER component of Q_{r} Q_{m} FOURIER component of Q_{} M_{rm} FOURIER component of M_{r} M_{m} FOURIER component of M_{} M_{rm} FOURIER component of M_{r} θ_{m} rotational angle of disk or cylindrical shell of pulley γ_{m} disk bending state variable vector of FOURIER component m m_{rm} disk circumferential harmonic resultant bending moment V_{rm} disk circumferential harmonic resultant transverse shear force A_{m} matrix containing coefficients of governing ODEs for disk bending deformation of FOURIER component m B_{m} vector representing the nonhomogeneous term of the governing ODEs for disk bending deformation of FOURIER component m U^{BD}_{m} TMB disk bending element displacement vector K_{BDextm} TMB disk bending element stiffness matrix F^{BD}_{extm} TMB disk bending element external force vector F^{BD}_{intm} TMB disk bending element internal force vector N_{r} normal force acting on the disk cross section perpendicular to radial direction N_{} normal force acting on the disk cross section perpendicular to circumferential direction N_{r} inplane shear force acting on the disk cross sections perpendicular to disk radial direction and disk circumferential direction N_{rm} FOURIER component of N_{r} N_{m} FOURIER component of N_{} N_{rm} FOURIER component of N_{r} η_{m} disk planestress state variable vector of FOURIER component m C_{m} matrix containing coefficients of governing ODEs for disk Inplane deformation of FOURIER component m U^{PN}_{m} TMB disk planestress element displacement vector K_{PNm} TMB disk planestress element stiffness matrix F^{PN}_{extm} TMB disk planestress element external force vector F^{PN}_{intm} TMB disk planestress element internal force vector U^{DK}_{m} TMB disk element displacement vector K_{DKm} TMB disk element stiffness matrix F^{DK}_{extm} TMB disk element external force vector F^{DK}_{intm} TMB disk element internal force vector R cylindrical shell radius N_{1} normal force acting on the cylindrical shell cross section perpendicular to pulley axial direction N_{2} normal force acting on the cylindrical shell cross section perpendicular to pulley circumferential direction S membrane shear force acting on the cylindrical shell cross sections perpendicular to pulley axial direction and circumferential direction Q_{1} transverse shear force acting on the cylindrical shell cross section perpendicular to pulley axial direction Q_{2} transverse shear force acting on the cylindrical shell cross section perpendicular to pulley circumferential direction M_{1} bending moment acting on the cylindrical shell cross section perpendicular to pulley axial direction M_{2} bending moment acting on the cylindrical shell cross section perpendicular to pulley circumferential direction M_{12} twisting moment acting on the cylindrical shell cross sections perpendicular to pulley axial direction and circumferential direction f_{z} external load acting on the neutral surface of shell in pulley axial direction f_{} external load acting on the neutral surface of shell in pulley circumferential direction f_{r} external load acting on the neutral surface of shell in pulley radial direction f_{zm} FOURIER component of f_{z} f_{m} FOURIER component of f_{} f_{rm} FOURIER component of f_{r} N_{1m} FOURIER component of N_{1} N_{2m} FOURIER component of N_{2} M_{1m} FOURIER component of S M_{2m} FOURIER component of M_{1} M_{12m} FOURIER component of M_{2} V_{1m} FOURIER component of M_{12} x_{m} FOURIER component of element boundary equivalent transverse shear force J_{m} cylindrical shell state variable vector of FOURIER component m I_{m} matrix containing coefficients of governing ODEs for cylindrical shell of FOURIER component m M_{12} vector representing the nonhomogeneous term of the governing ODEs for cylindrical shell of FOURIER component m U^{CS}_{m} TMB cylindrical shell element displacement vector K_{cs m} TMB cylindrical shell stiffness matrix F^{CS}_{extm} TMB cylindrical shell element external force vector F^{CS}_{intm} TMB cylindrical shell element internal force vector Ψ_{m} state variable vector of FOURIER component m _{S} axial or radial coordinate H_{m} matrix containing coefficients of governing ODEs of a pulley component L_{m} vector representing the nonhomogeneous term of the governing ODEs of a pulley component U_{m} generalized displacements of state variable vector of a pulley component F_{m} generalized forces of state variable vector of a pulley component T_{m} transfer matrix for the governing ODEs of a pulley component D_{m} TMB element displacement vector of a pulley component K_{m} TMB element stiffness matrix of a pulley component F_{int m} TMB element internal force vector of a pulley component F_{ext m} TMB element external force vector of a pulley component 1. Introduction
An engineered class belt conveyor pulley typically consists of a cylindrical shell, two end disks with variable thickness, a shaft, and two locking devices connecting end disks to the shaft as shown in Fig, 1. The pulley is usually subjected to severe bending due to very high belt tensions and locking assembly pressures, In the design of such a pulley, it is necessary to take into account the possibility of fatigue failure. Costly failures in large conveyor pulleys have led designers to seek detailed stress fatigue or endurance analysis. To date, two types of approaches for pulley stress analysis have been reported in the literature. One is the classical mechanics approach developed by LANGE [1] and SCHMOLTZI [2]. The other is the finite element method (FEM) employed by VODSTRICL [3], DANIEL [4] and SETHI et al. [5]. Both types of approaches have advantages and disadvantages.
Fig. 1: Crosssection of pulley assembly
The classical mechanics approach developed by LANGE and SCHMOLTZI is an approximate analytical approach, providing a closedform solution for stresses in a pulley. The advantages of this method are that it is easy to program and takes a very short execution time to obtain a solution, The disadvantage is that the stress solution is not accurate at the locations near the connection region between the shell and end disks because of its poor approximation in treating the elastic coupling between these components. Specifically, the displacement of the end disk and shell are not coupled at their connection. This leads to significant errors in the stress and strain field about the connectors as will be shown.
The FEM has just the opposite advantages and disadvantages of the LANGE classical method. The major advantage of FEM is its ease of treating complex geometry and boundary conditions. The major disadvantage is its long execution time coupled with its need for an experienced user to generate a proper finite element mesh.
In this paper, a new method called the Modified Transfer Matrix (MTM) method is presented, This method circumvents the disadvantages of both the LANGE classical method and the conventional FEM. The MTM method proves to be a very effective and efficient approach in providing an accurate pulley assembly stress solution for any loading condition pulley.
Historically, the transfer matrix method was developed several decades ago [6] and was very popular in solving one dimensional static and dynamic problems before the advent of the FEM. Even today, this method is still useful in providing dosed form solutions to certain elasticity problems with simple boundary conditions [7]. Although there is a limitation in handling complicated boundary conditions, such as the boundary conditions for a pulley, the solution obtained by using the transfer matrix method is exact. In this paper, it is shown that the limitation of the transfer matrix method can be overcome if the transfer matrix is reformulated by using finite element concepts. The reformulated transfer matrix is essentially a special finite element, The new method using these special finite elements, called transfer matrix based (TMB) finite elements, is capable of solving a class of structural elasticity problems (including the elasticity problem of a pulley), whose governing differential equations can be reduced to a set of ordinary differential equations (ODEs). Regardless of how few of these TMB finite elements are used in a model, the solution obtained by this MTM method is generally very accurate due to the nature of the transfer matrix method.
Based on the MTM method, a computer program for pulley stress analysis named PSTRESS 3,0 has been developed, This program can provide stress solutions and perform fatigue analysis for most pulleys, with the characteristic geometry shown in Fig. 1, The pulley can be subjected to any type of nonuniform surface pressure, sheer loading, and prescribed locking pressure.
In section 2, the general ideas for deriving TMB elements for beam (i.e. shaft), enddisk plate with variable thickness, and cylindrical shell are presented. In section 3, the assembly of these elements to model a pulley in PSTRESS 3.0 is discussed. In section 4, an example of a belt conveyor pulley is numerically solved by PSTRESS 3.0 and the results are compared with those obtained using a finely meshed FEM (ANSYS) solution.
2. TMB Finite Elements for Beam, Disk Plate and Cylindrical Shell
Stresses and displacements in a pulley can be expressed in terms of FOURIER series with respect to the circumferential angle because of the pulley's axisymmetric geometry. Each FOURIER component of the solutions can be determined by solving a set of corresponding governing differential equations, which are uncoupled with the governing equations for other Fourier components. In Appendices A, B and C, it is shown that the governing equations for Fourier components for shaft, enddisk plate with nonuniform thickness, arid cylindrical shell of a pulley can be reduced to a set of ODEs of the first order respectively. The general solutions to these ODEs can be expressed in terms of the transfer matrix. By following the procedure described in Appendix D, the general solutions can be reorganized in a finite element form as below
K_{m}D_{m} = F_{int m} + F_{ext m} (1)
where K_{m} is the TMB element stiffness matrix, D_{m} is the element displacement vector, F_{int m} is the element internal force vector, F_{ext m} is the element external force vector, and the subscript "m" denotes the FOURIER component number, The detailed procedures of developing TMB elements for shaft, enddisk and cylindrical shell are given in Appendices A, B and C, respectively.
Remarks;
As seen in the above discussion, the general solution to the governing differential equations for a pulley can be finally transformed into a finiteelement form. This allows us to exploit many finite element analysis (FEA) capabilities to resolve pulley stresses using our MTM method, The most valuable FEA capability to be employed is the way of treating complicated boundary conditions. Therefore, using the MINI method, we can easily take into account the elastic coupling between the rim and the enddisk by following FEA assembly procedures, and implement the locking assembly pressure by using the FEA approach of treating mechanical Interference between two bodies. Both of these problems cannot be easily or precisely handled by most classical methods.
3. Assembly of TMB Elements for a Pulley Model in PSTRESS 3.0
In PSTRESS 3.0, the above derived TMB beam, disk plate and cylindrical shell element stiffness matrices are brought together to font a global stiffness matrix for a pulley in essentially the same way as that in conventional FEM, However, care must be taken at two locations, where special element assembly methods are required.
The first location is the connection region between the shell and the disk shown in Fig, 2a, where the finite dimension of the joint has to be taken into account in the stiffness matrix. The conventional way of treating the whole region as a single node would cause significant error in stress solutions near this region. One of the reasonable ways of providing correct elastic stiffness to connect the rim and the disk is treating the whole region as a special element by applying a substructure method. In PSTRESS 3.0, such a special element shown in Fig. 2b is developed.
Fig. 2: Connecting element between rim and enddisk
The second location is the connection point between the locking device and the shaft, One thing that must be kept in mind when assembling elements in this location is that the final solution consists of many FOURIER components, of which the shaft element only contributes to the FOURIER components of m = 0,1 and 1. In fact, in a pulley structure, the shaft deformation is governed only by these three FOURIER components due to its slender geometry. For the component of m = 0, the corresponding shaft deformation can be exactly modeled by a TMB cylindrical shell element presented in Appendix C. When a TMB beam element described in Appendix A is used, representing the shaft deformation of components of m = 1 and 1 (i.e., bending deformation), the following constraints on the deformation of the connecting point must be imposed:
W_{m} = m V_{m}  (2) 
U_{m} = r_{shaft}θ_{m}  (3) 
W_{s} = W_{m}  (4) 
θ_{s} = θ_{m}  (5) 
W_{s} and θ_{s} are the shaft deflection and rotational angle respectively at the connecting point, where the subscript s denotes shaft deformation. The general definitions of W_{s} and θ_{s} are given in Appendix A.
U_{m}, V_{m}, W_{m} and θ_{m} are disk plate displacements and rotational angle at the connecting point, where the subscript m denotes the FOURIER component number: m = 1 or 1. The general definitions of U_{m}, V_{m}, W_{m} and θ_{m}, are given in Appendix B. r_{shaft} is the shaft radius at The connection point.
Such constraints are easy to implement in a finite element model by using FEA static condensation or penalty methods [10]. In PSTRESS 3.0, the static condensation method is employed. Finally, it must be pointed out that the TMB beam element stiffness matrix and corresponding nodal forces must be multiplied by a factor of 2 before they are assembled into global equations, due to the difference between actual forces and harmonic forces.
4. Numerical Example
Consider the belt conveyor pulley shown in Fig, 1, which is supported on two bearings and subjected to a locking pressure of 115.71 mPa (16,778 psi) at the interface between the locking device and disk hub. The shell circumferential surface pressure and shear loading between circumferential angles of 83 and 254% are developed from unequal belt tensions T_{1} = 1,017.8 kN (228,800 lb) and T_{2} = 632.98 kN (142.300 lb) shown in Fig. 3. The material properties and geometrical parameters of the pulley are given in Table 1, This pulley is analyzed by using PSTRESS 3.0, ANSYS 4.4 and CDIs derivation of LANGES classical method, respectively. Because of symmetry, only one quarter of pulley cross section is modeled.
Fig. 3: Load on pulley due to belt tension
In the PSTRESS model, the rim is modeled with 2 TMB elements, the disk is modeled with 6 TMB elements, the shaft is modeled with 3 TMB elements, and 71 FOURIER components are used. The reason for using more TMB elements for the disk is the necessity of taking account of the nonuniform thickness of the disk. In the ANSYS FEA model, 5,000 axisymmetric structural solid elements (with non axissymmetric loading) are employed in the 20 crosssection. The use of the ANSYS FEM package to analyze a pulley is discussed in [5].
Material Property  
Young's modulus, MPSI  30  
Poisson's ratio  0.3  
Rim Geometry  
Rim length, inches  82  
Rim outer diameter, inches  54  
Rim thickness, inches  1.5  
Belt width, inches  72  
Disk Geometry  
Locking device width, inches  3.3  
Hub outer diameter, inches  27.6  
Hub inner diameter, inches  20.27  
Hub width, inches  6.7  
Fillet radius at hub, inches  3.67  
Filet radius at rim, inches  1.20  
Disk thickness between hub and rim, inches  


Shaft Geometry  
Diameter, inches  16.535  
Shaft length, inches  121  
Distance between bearing centers, inches  104  
Distance between disk centers, inches  74.13 
Table 1: Material properties and geometrical parameters.
Figs. 47 show the PSTRESS numerical results compared with the ANSYS results. From these figures, it is seen that at location A of the rim and location D of the disk the agreement between the results Of the MTM method and the results of the conventional FEM is good. At location B of the rim and location C of the disk the agreement is still good, but some inaccuracy is observed, The reason may be that locations B and C are within the connection region between the rim and disk, where the 3D stress state is more significant and cannot be fully taken into account in 2D shell and plate theories. According to St. VENANTS principle and our experience, This 3D stress state has only a very localized effect on pulley stress solution when the thicknesses of rim and disk are relatively small compared with the length and radius of the rim. It must be noted that the MTM method is much more efficient than the conventional FEM. PSTRESS 3.0 takes approximately 30 seconds to obtain a solution on an IBM PC 486, including the fatigue analysis. The FEM (ANSYS) solution takes 1224 hours on an IBM RISC 6000 workstation.
Fig. 4: Stresses at Location A (inside of rim) PSTRESS 3.0 Analysis 
Fig. 5: Stresses at Location B (inside of rim) PSTRESS 3.0 Analysis 
Fig. 6: Stresses at Location C (inside of disk) PSTRESS 3.0 Analysis 
Fig. 7: Stresses at Location D (inside of disk) PSTRESS 3.0 Analysis 
Figs. 811 show the comparison of numerical results between LANGES solution and the FEM solution. Except at location A, LANGE's solution does not agree with the FEM solution. The poor agreement is due to the errors in treating the elastic coupling between the rim and the enddisk in LANGEs method.
Fig. 8: Stresses at Location A (Inside of rim) LANGE Analysis 
Fig. 9: Stresses at Location B (Inside of rim) LANGE Analysis 
Fig. 10: Stresses at Location C (Inside of disk) LANGE Analysis 
Fig. 11: Stresses at Location D (Inside of disk) LANGE Analysis 
Figs. 12 and 13 show both the PSTRESS and ANSYS results at two corners of the interface between the locking device and the shaft (locations EE of Fig. 1).At these two locations, 3D stress state is much more significant than at locations B and C. In order to produce more accurate stress solutions at the two corners, we introduce stress concentration factors in zero and first order FOURIER component solutions by using our empirical formulae built in PSTRESS 3.0. As seen in Figs. 12 and 13, these corrected solutions agree well with ANSYS solutions.
Fig. 12: Stresses at Location E (inside of pulley) PSTRESS 3.0 Analysis 
Fig. 13: Stresses at Location E (outside of pulley) PSTRESS 3.0 Analysis 
5. Conclusions
A new pulley stress analysis method which is based on reformulated transfer matrix has been developed. An accurate solution can be obtained by using this method. Three transfermatrixbased elements for the shaft, disk plate and cylindrical shell have been developed. A numerical example has been given, which demonstrates the merits of this new method.
References
 LANGE, H.: Investigations on Stress in Belt Conveyor Pul leys; Doctoral thesis, Technical University Hannover, 1963.
 SCHMOLTZI, W.: The Design of Conveyor Belt Pulleys with Continuous Shafts; Doctoral thesis, Technical University, Hannover 1974.
 V0DSTRCIL, P.: Analysis of Belt Conveyor Pulley Using Finite Element Method; Proc. 4th i nt. Conf. in Australia on Finite Element Methods, University of Melbourne, Aug. 1820, 1982.
 DANIEL, W.J.T.: Development of a Conveyor Pulley Stress Analysis Package; Proc. nt. Conf. on Bulk Material Storage, Handling and Transportation, Newcastle, Aug. 2224, 1983.
 SETHI, V. and NORDELL, L.K.: Modern Pulley Design Techniques and Failure Analysis Methods; Proceedings of SME Annual Meeting & Exhibit, Peno Nevada, USA, Feb. 1518,1993.
 PESTEL, S.C. and LECKIE, F.A.: Matrix Methods in Elastomechanics; New York McGrawHill Book Co., 1963.
 YEH KAIYUAN: General Solution on Certain Problems of Elasticity with NonHomogeneity and Variable Thickness, The Advances of Applied Mathematics and Mechanics, Vol. 1; China Academic Publishers, pp 240273, 1987.
 BOYCE, W.E. and DIPRIMA, R.C.: Elementary Differential Equations; Fifth Edition, John Wiley & Sons Inc., New York, 1992.
 TIMOSHENKO, S. and WOINOWSKYKRIEGER. S.: Theory of Plates and Shells; McGrawHill Book Co., Second Edition, 1959.
 COOK, R.D., MALKUS, D.S. and PLE5HA, M.E.: Concepts and Applications of Finite Element Analysis; Third Edition, John Wiley & Sons Inc., 1969.
Appendix A
TMB Finite Element for Beam
Fig. A1 shows forces that act on a differential beam. Loads P, M, and q are shown in their positive sense. z is the axial coordinate.
Fig. A1: Forces that act on a differential element of beam
The equilibrium equations are
dP/dz = q (A1)
dM/dz = P (A2)By using TIM0SHENKOs beam theory, we have
dW_{s}/dz = θ_{s} + P/kaG (A3)
dθ_{s}/dz = M/EI (A4)where the subscript s denotes the shaft deformation, W_{s} is the shaft neutral axis displacement, θ_{s} is the slope due to bending, dW_{s}/dz is the slope of the center line of the beam, k is a shape factor equal to 0.75 for circular cross section, EI is the bending stiffness, a is the crosssection area, and G is the shear modulus.
Eqs. (A1)(A4) can be written in a matrix form
dβ/dz = Aβ + B (A5)
where
β = (W_{s}, θ_{s}, P, M)^{T} (A6)
B = (0, 0, q, 0)^{T} (A7)and
A = 
[ 
0  1  1/kaG  0  ]  (A8) 
0  0  0  1/EI  
0  0  0  0  
0  0  1  0 
where β is called state variable vector.
According to the ODE theory [8], the general solution to Eq. (A5) can be expressed as
β(z) = T(z)β(z_{0}) + T(z)  ^{z} ∫ _{z0} 
T(s)^{1}B(s)ds (A9) 
where z_{0} ≤ z ≤ z_{L}, z_{0} and z_{L} are the coordinates corresponding to two ends of the beam, and T(z) is the transfer matrix satisfying
dT/dz = AT and T(z_{0}) = l (A10)
where l is the identity matrix. Following the general procedure described in [8]. we can obtain the following closedform transfer matrix
T(z) =  [  1  zz_{0}  zz_{0}  +  (zz_{0})  (zz_{0})  ]  (A11) 
kaG  2EI  2EI  
0  1  (zz_{0})  zz_{0}  
2EI  EI  
0  0  1  0  
0  0  zz_{0}  1 
where z_{0} ≤ z ≤ z_{L}. Following the procedure described in Appendix D, we can obtain the following finite element equation, which is equivalent to Eq. (A9)
K_{BM}U_{BM} = F^{BM}_{int} + F^{BM}_{ext} (A12)
where
U_{BM} = (W_{s} (Z_{L}), θ_{s}(Z_{L}), W_{s}(Z_{0}), θ_{s}(Z_{0}))^{T} (A13)
F^{BM}_{ext} is the beam element external force vector, F^{BM}_{int} is the beam element internal force vector, and K_{BM} is the beam element stiffness matrix, which can be expressed:
KBM =  [  12EI  Symmetric  ]  (A14)  
L(1+α)  
6EI  (4+α)EI  
L(1+α)  L(1+α)  
12EI  6EI  12EI  
L(1+α)  L(1+α)  L(1+α)  
6EI  (2α)EI  6EI  (a+α)EI  
L(1+α)  L(1+α)  L(1+α)  L(1+α) 
Where
α = 12EI/LkaG and L = Z_{L}  Z_{0} (A15)
Remarks:
It is not surprising to note that the TMB finite element stiffness matrix for a beam derived as Eq. (A14) is identical to the beam stiffness matrix derived by using conventional finite element method because the polynomialtype trial function used in conventional FEM exactly represents the actual beam displacement. However, for other types of structures such as plate and shell, the TMB element stiffness matrices may not be the same as the conventional finite element stiffness matrices.
Appendix B
TMB Finite Element for Disk Plate with Variable Thickness
B.1 TMB Bending Element
Fig. B1 shows a differential element of plate subjected to bending loads. x, y and z are Cartesian coordinates with z coinciding with pulley axial direction, x horizontal direction, and y vertical direction. r and are disk cylindrical coordinates. Loads Q_{r}, Q_{}, M_{r}, M_{} and M_{r} are shown in their positive sense.
Fig. B1: Bending forces that act on a different element of plate
The equilibrium equations are:
Q_{r} =  ∂M_{r}  +  M_{r}M_{}    1  ∂M_{r}  (B1)  
∂r  r  r  ∂  
Q_{} =  ∂M_{r}  +2  M_{r}    1  ∂M_{}  (B2)  
∂r  r  r  ∂  
∂Q_{r}  +  Q_{r}    1  ∂Q_{}  = g(r, )  (B3)  
∂r  r  r  ∂ 
where g is the external transverse force acting on the neutral surface of the plate, and equations representing the forcedeformation relationship are:
M_{r} = D  [  ∂u  +  (  1  ∂u  +  1  ∂u  )]  (B4) 
∂r  r  ∂r  r  ∂  
M_{} = D  [  ∂u  +  1  ∂u  +  1  ∂u  ]  (B5)  
∂r  r  ∂r  r  ∂  
M_{r} = D(1) 
[  1  ∂u    1  ∂u  ]  (B6)  
r  ∂r∂  r  ∂ 
where U is the plate neutral surface transverse displacement (see Fig. B1), is the POISSON's ratio,
D = Et/12(1  ) (B7)
E is the YOUNG's modulus, t is the thickness of the plate, which is a function of r, Within the element, we assume
t = cr^{p} (B8)
where a and p are constants. Also, we assume the following FOURIER series for the components of displacement and forces
u =  ∞ Σ m=0 
u_{m}cosm +  ∞ Σ m=1 
u_{m}sinm (B9) 
g =  ∞ Σ m=0 
g_{m}cosm +  ∞ Σ m=1 
g_{m}sinm (B10) 
Q_{r} =  ∞ Σ m=0 
Q_{rm}cosm +  ∞ Σ m=1 
Q_{r,m}sinm (B11) 
Q_{} =  ∞ Σ m=1 
Q_{m}sinm +  ∞ Σ m=0 
Q_{,m}cosm (B12) 
M_{r} =  ∞ Σ m=0 
M_{rm}cosm +  ∞ Σ m=1 
M_{r,m}sinm (B13) 
M_{} =  ∞ Σ m=0 
M_{m}cosm +  ∞ Σ m=1 
M_{,m}sinm (B14) 
M_{r} =  ∞ Σ m=1 
M_{r}sinm +  ∞ Σ m=0 
M_{r,m}cosm (B15) 
where u_{m}, P_{m}, Q_{rm}, Q_{m}, M_{rm}, M_{m}, M_{rm }(m=O.1,2,...) are functions of r only. Substituting Eqs. (89)(B15) into Eqs. (B1)(B6), we have the following ordinary differential equations
Q_{rm} = M'_{rm} +  1  (M_{rm}  M_{m})   m  M_{rm}  (B16)  
r  r  
Q_{m} = M'_{rm} +  2  M_{rm} +  m 
M_{m} 
(B17)  
r  r  
Q'_{rm} + 
1  Q_{rm}   m  Q_{m} = g_{m}  (B18)  
r  r  
M_{rm} = D  [  u"_{m} +  (  1  u'_{m}   m  u_{m}  )]  (B19) 
r  r  
M_{m} =  D  [  u"_{m} +  1  u'_{m}   m  u_{m} 
] 
(B20)  
r  r  
M_{rm} = D(1) 
[ 
  m  u'_{m} +  m  u_{m} 
] 
(B21)  
r  r 
where m = 0, 1, 2, ..., and the prime represents derivative with respect to r. Introducing the following state variables
γm = (u_{m}, θ_{m}, r_{m}, M_{rm})^{T} (B22)
Where
θ_{m} = u'_{m}  (B23) 
m_{rm} = 2πr M_{rm}  (B24) 
V_{rm} = 2πr (Q_{m}  m/r M_{rm})  (B25) 
Eliminating five variables, M_{rm}, M_{m}, M_{rm}, Q_{rm} and Q_{m} among Eqs. (B16)(B21) and (B23)(B25), we can obtain the following matrixform equation
dγ_{m}/dr = A_{m} γ_{m} + B_{m} (m = 0, 1, 2, ...) (B26)
Where Am =
[  0  1  0  0  ]  (B27) 
m  _  0  1  
r  r  2πDr  
2πDm(2  2 + m  m)  2πDm(3  2  )  _ 1  m  
r  r  r  r  
2πDm(3 + 2 + )  2πD(1  + 2m  2m)  1  (1  )  
r  r  r 
B_{m} = (0, 0, 2πrg_{m}, 0)^{T} (B28)
Therefore, following the procedure described in Appendix D, we can obtain the following finite element equation for a circular plate ring (r_{0} ≤ r ≤ r_{L}) with variable thickness and subjected to harmonic bending load
K_{BD m}U^{BD}_{m} = F^{BD}_{int m} + F^{BD}_{ext m} (B29)
Where
U^{BD}_{m} = (U_{Lm}, _{Lm}, U_{om}, θ_{0m})^{T} (B30)
K_{BD m} is the plate element bending stiffness matrix, F^{BD}_{ext m} is the element external force vector, F^{BD}_{int m} is the element internal force vector, the subscript L and 0 denote locations at r = r_{L} and r_{0}, respectively, and m denotes the FOURIER component number. Due to the complexity in deriving the transfer matrix for Eq. (B26), it is much more difficult to obtain a closedform expression for K_{BD m} of Eq. (B29) than for KBM of Eq. (A12). Instead, we can very accurately calculate K_{BD m} by using our computer program of PSTRESS 3.0 mentioned in Remark ii of Appendix D.
B.2 TBM Plane Stress Element
Fig. B2 shows a differential element of plate subjected to inplane loading. Loads N_{r}, N_{}, and N_{r} are shown in their positive sense.
Fig. B2: Inplane forces that act on a differential element of plate
The definitions of x, y, Z, r and are the same as in section B.1. The equilibrium equations are
∂N_{r}  +  N_{r}  N_{}  +  1  ∂N_{r}  = 0 (B31) 
∂r  r  r  ∂  
∂N_{r}  + 2  N_{r}  +  1  ∂N_{}  = 0 (B32) 
∂r  r  r  ∂ 
The equations for forcedeformation relationship are
N_{r} =  Et  (  ∂w  +  1  ∂v  +  w  )  (B33)  
(1  )  ∂r  r  ∂  r  
N_{} =  Et  (  ∂w  +  1  ∂v  +  w  )  (B34)  
(1  )  ∂r  r  ∂  r  
N_{r} =  Et  (  ∂v  +  1  ∂w    v  )  (B35)  
2(1 + )  ∂r  r  ∂θ  r 
where the definitions of , E and t are the same as in section B.1, arid v and w are the displacements of the neutral surface in circumferential and radial directions, respectively.
We assume the following FOURIER series for the components of displacements and forces
v =  ∞ Σ m=1 
v_{m}sinm +  ∞ Σ m=0 
v_{m}cosm (B36) 
w =  ∞ Σ m=0 
w_{m}cosm +  ∞ Σ m=1 
w_{m}sinm (B37) 
N_{r} =  ∞ Σ m=0 
N_{rm}cosm +  ∞ Σ m=1 
N_{r,m}sinm (B38) 
N_{} =  ∞ Σ m=0 
N_{m}cosm +  ∞ Σ m=1 
N_{,m}sinm (B39) 
N_{r} =  ∞ Σ m=1 
N_{rm}sinm +  ∞ Σ m=0 
N_{rm}cosm (B40) 
where v_{m}w_{m}, N_{rm}, N_{m} and N_{rm}(m = 0,1,2, ...) are functions of r only. Substituting Eqs. (B36)(B40) into Eqs. (B31)(B35), we have the following ordinary differential equations
N'_{rm} +  1  (N_{rm}  N_{rm})+  m  N_{rm} = 0  (B41)  
r  r  
N'_{rm} +  2 
N_{rm}  
m  N_{m} = 0  (B42)  
r  r  
N_{rm} =  Et  [  w'_{m} +  (  m  v_{m} +  1  wm  )]  (B43) 
(1  )  r  r  
N_{m} =  Et  [  w'_{m} +  m  v_{m} +  1  wm  ]  (B44)  
(1  )  r  r  
N_{rm} =  Et  [  v'_{m}   m  w_{m}   1  vm  ]  (B45)  
2(1 + )  r  r 
where m = 0, 1, 2, 3,.... Introducing the following state variables
η_{m} = (w_{m}, v_{m}, 2πrN_{rm} , 2πrN_{rm})^{T} (B46)
and eliminating N_{m} among Eqs. (B41)(B45). we obtain
dη_{m}/dr = C_{m} η_{m} (m = 0, 1, 2, ...) (B47)
where
C_{m} =  [  _  _ m  1  0  ]  (B48) 
r  r  2πEtr  
m  1  0  1 +  
r  r  2πEtr  
2πEt  2πEtm  _ 1   _ m  
r  r  r  r  
2πEtm  2πEtm  m  _ 2  
r  r  r  r 
Thus, following the similar procedure in deriving Eq. (B29) , we can finally derive the following TMB finite element equation for a circular ring (r_{0} ≤ r ≤ r_{L}) with variable thickness, subjected to inplane harmonic loading
K_{PN m} U^{PN}_{m} = F^{PN}_{int m} + F^{PN}_{ext m} (B49)
where
U^{PN}_{m} = (W_{Lm}, V_{Lm}, W_{0m}, V_{0m})^{T} (B50)
K_{PN m} is the element planestress stiffness matrix, F^{PN}_{ext m} is the element external force vector, F^{PN}_{int m} is the element internal force vector, and the subscript L and C denote locations at r = r_{L} and r_{0} respectively, and m denotes the FOURIER component number.
Combining Eq. (B29) and Eq. (B49), we can form a TMB shell element for disk plate, which is subjected to both bending and inplane loading
K_{DK m} U^{DK}_{m} = F^{Dk}_{int m} + F^{DK}_{ext m} (B51)
where K_{DK m} is the disk element stiffness matrix, U^{DK}_{m} is the element displacement vector, F^{DK}_{ext m} is the element extemal force vector, and F^{Dk}_{int m} is the element internal force vector. They can be calculated by the following formulae
U^{DK}_{m} = (U_{Lm}, V_{Lm}, W_{Lm}, θ_{Lm}, U_{0m}, V_{0m}, W_{0m}, θ_{0m})^{T}  (B52) 
K_{DK m} = S^{T}_{1}K_{BD m}S_{1 }+ S^{T}_{2}K_{PN m}S_{2}  (B53) 
F^{DK}_{int m} = S^{T}_{1}F^{BD}_{int m }+ S^{T}_{2}F^{PN}_{int m}  (B54) 
F^{DK}_{ext m} = S^{T}_{1} F^{BD}_{ext m } + S^{T}_{2} F^{PN}_{ext m}  (B55) 
S_{1} =  [  1  0  0  0  0  0  0  0  ]  (B56) 
0  0  0  1  0  0  0  0  
0  0  0  0  1  0  0  0  
0  0  0  0  0  0  0  1 
Remarks:
 Any other state variables than those defined in Eq. (B22) and Eq. (B46) cannot be employed because they may lead to incorrect element stiffness matrices and force vectors due to the violation of MAXWELL'S reciprocal theorem.
 Stiffness matrices obtained from transfer matrices must be symmetric. Any mistakes in choosing state variables, in deriving equations, and/or in numerical programming may lead to asymmetric stiffness matrices.
 If the stress calculation is desired, we can follow the following steps
 calculate values of state variables at element nodes,
 use transfer matrix and Eq. (D3) to calculate values of state variables at any interior location of the element, where stress values are desired,
 calculate first derivatives of state variables at desired locations by using Eq. (B26) and Eq. (B47),
 calculate M_{rm}, M_{m}, M_{r}_{}_{m}, N_{rm}, N_{m} and N_{rm} by using Eqs. (B19)(B21) and Eqs. (B43)(B45), and
 calculate FOURIER components of stresses by using calculated forces in step d.
Appendix C
TMB Finite Element for Cylindrical Shell
Fig. C1 shows forces acting on a differential element of a cylindrical shell with constant thickness t and radius R, z and are cylindrical coordinates. Loads N_{1}, N_{2}, S, Q_{1}, Q_{2}, M_{1}, M_{2} and M_{12} are shown in their positive sense.
Fig. C1: Forces that act on a differential element of cylindrical
The equilibrium equations, according to [9] are
∂N_{1}  +  ∂S  + f_{z} = 0  (C1)  
∂z  R∂  
∂N_{2}  +  ∂S  + f_{} = 0  (C2)  
R∂  ∂z  
_ N_{2}  +  ∂Q_{1}  +  ∂Q_{2}  + f_{r} = 0  (C3) 
R  ∂z  R∂  
Q_{2} =  ∂M_{12}  +  ∂M_{2}  (C4)  
∂z  R∂  
Q_{1} =  ∂M_{1}  +  ∂M_{12}  (C5)  
∂z  R∂ 
where f_{z}, f_{} and f_{r} are external loads acting on the neutral surface of the shell in axial, circumferential and radial directions respectively, and the equations for forcedisplacement relationship, according to [9]. are
N_{1} =  Et  (  ∂u  +  ∂v  +  w  )  (C6)  
(1)  ∂z  R∂  R  
N_{2} =  Et  (  ∂u  +  ∂v  +  w  )  (C7)  
(1)  ∂z  R∂  R  
S =  Et  (  ∂u  +  ∂v  )  (C8)  
(12)  R∂  ∂z  
M_{1} = D  (  ∂w  +  ∂u  )  (C9)  
∂z  R∂  
M_{2} = D  (  ∂w  +  ∂u  )  (C10)  
∂z  R∂  
M_{12} = (1  )D  ∂w  (C11)  
R∂z∂ 
where u, v and w are shell neutral surface displacements in axial, circumferential and radial directions, respectively,
D = Et/12(1  )
the definitions of E and are the same as in section B, and t is the constant shell thickness. We assume the following FOURIER components of displacements and forces
f_{z} =  ∞ Σ m=0 
f_{zm}cosm +  ∞ Σ m=1 
f_{z,m}sinm  (C12) 
f_{} =  ∞ Σ m=1 
f_{m}sinm +  ∞ Σ m=0 
f_{,m}cosm  (C13) 
f_{r} =  ∞ Σ m=0 
f_{rm}cosm +  ∞ Σ m=1 
f_{r,m}sinm  (C14) 
u =  ∞ Σ m=0 
u_{m}cosm +  ∞ Σ m=1 
u_{m}sinm  (C15) 
v =  ∞ Σ m=1 
v_{m}sinm +  ∞ Σ m=0 
v_{m}cosm  (C16) 
w =  ∞ Σ m=0 
w_{m}cosm +  ∞ Σ m=1 
w_{m}sinm  (C17) 
N_{1} =  ∞ Σ m=0 
N_{1m}cosm +  ∞ Σ m=1 
N_{1,m}sinm  (C18) 
N_{2} =  ∞ Σ m=0 
N_{2m}cosm +  ∞ Σ m=1 
N_{2,m}sinm  (C19) 
S =  ∞ Σ m=1 
S_{m}sinm +  ∞ Σ m=0 
S_{m}cosm  (C20) 
M_{1} =  ∞ Σ m=0 
M_{1m}cosm +  ∞ Σ m=1 
M_{1,m}sinm  (C21) 
M_{2} =  ∞ Σ m=0 
M_{2m}cosm +  ∞ Σ m=1 
M_{2,m}sinm  (C22) 
M_{12} =  ∞ Σ m=1 
M_{12m}sinm +  ∞ Σ m=0 
M_{12,m}cosm  (C23) 
Q_{1} =  ∞ Σ m=0 
Q_{1m}cosm +  ∞ Σ m=1 
Q_{1,m}sinm  (C24) 
Q_{2} =  ∞ Σ m=1 
Q_{2m}sinm +  ∞ Σ m=0 
Q_{2,m}cosm  (C25) 
where all FOURIER coefficients are functions of z only. Substituting Eqs. (C12)(C25) into Eqs. (C1)(C11), introducing the following state variables
ξ_{m} = (U_{m}, V_{m}, W_{m}, θ_{m}, 2πRN_{1m}, 2πRS_{m}, 2πRV_{1m}, 2πRM_{1m})^{T} (C26)
Where
θ_{m} = W'_{m} (C27)
and
V_{1m} = Q_{1m} + m/R M_{12m} (C28)
and eliminating five variables, N_{2m}, M_{2m}, M_{12m}, Q_{1m} and Q_{2m }we can finally obtain the following matrixform equation after a long tedious procedure of mathematical derivation
dξ_{m}/dz = J_{m}ξ_{m} + I_{m} (m = 0, 1, 2, ...) (C29)
Where
J_{m} =  [  0  _ m  _  0  1   0  0  0  ]  (C30)  
R  R  2πEtR  
m  0  0  0  0  1 +  0  0  
R  πEtR  
0  0  0  1  0  0  0  0  
0  0  _ m  0  0  0  0  6(1  )  
R  πEtR  
0  0  0  0  0  _ m  0  0  
R  
0  πEtm  2πEtm  0  m  0  0  0  
R  R  R  
0  2πEtm  2πEt  (  1 +  tm^{4}  )  0  0  0  m  
R  R  12R  R  R  
0  0  0  πEtR  0  0  1  0  
3(1+)R 
and
I_{m} = (0, 0, 0, 0,  2πRf_{zm}, 2πRf_{m}, 2πRf_{rm}, 0)^{T} (C31)
Therefore, according to the conclusions of Appendix D, the transfer matrix for Eq. (C29) can be obtained by the computer program developed in PSTRESS 3.0, Based on the transfer matrix, the TMB finite element for cylindrical shell with z_{0} ≤ z ≤ z_{L} can be derived by following the procedure described in Appendix D, and the result can be expressed by
K_{CS m} U^{CS}_{m} = F^{CS}_{int m} + F^{CS}_{ext m} (C32)
where
U^{CS}_{m} = (U_{Lm}, V_{Lm}, W_{Lm}, θ_{Lm}, U_{0m}, V_{0m}, W_{0m}, θ_{0m})^{T} (C33)
K_{CS m} is the cylindricalshell stiffness matrix, F^{CS}_{ext m} is the element external force vector, F^{CS}_{int m} is the element internal force vector, and the subscripts L and 0 denote locations at z  z_{L}, and z = z_{0}, respectively, and m denotes the FOURIER component number.
If stress calculation for cylindrical shell is desired, we can follow the similar procedure described in Remark vi of section B.
Appendix D
Procedure of Deriving TMB Finite Element
In Appendices A, B and C, it is shown that the governing equations for FOURIER components for shaft, enddisk plate with nonuniform thickness, and cylindrical shell of a pulley can be reduced to a set of ODEs of first order respectively, and these ODEs can be finally unified as the following form
dψ_{m}/ds = H_{m}(s)ψ_{m} + L_{m} (m = 0, 1, 2, ...) (D1)
where the subscript m denotes the ordinary number of FOURIER components, H_{m} is a 2n x2n matrix. L_{m} is a 2n x 1 vector containing external loads, ψ_{m} is the state variable vector defined as
ψ_{m} = (U^{T}_{m}, F^{T}_{m})^{T} (D2)
U_{m} is an n x 1 vector containing generalized displacements, F_{m} an n x 1 vector containing corresponding generalized internal forces, and s is an axial or radial coordinate of the pulley, which is assumed within the range of s_{0} ≤ s ≤ s_{L}, where s_{0} and s_{L} are corresponding to two ends of a pulley component section (i,e, element).
Because Eq. (D1) is linear, the general solution is obtainable. By using the theory of ODE [8], the general solution to Eq. (D1) can be written as
ψ_{m}(S) = T_{m}(S) ψ_{m}(s_{0}) + T_{m}(S)  s ∫ s_{0} 
T_{m}(X)^{1}L_{m} (x)dx (D3) 
where T_{m}(s) is a 2n x 2n matrix called transfer matrix satisfying
dT_{m}/ds = H_{m}T_{m} and T_{m}(s_{0}) = I (D4)
where I is the identity matrix, For any set of linear ordinary differential equations, which can be expressed in the form of Eq. (D1), there exists a unique 2n x 2n transfer matrix T_{m}(s) satisfying Eq. (D4). If Eq. (D1) is a set of constant ODEs (i.e., H_{m} is independent of s) or it can be transformed into a set of constant ODEs, a closedform transfer matrix T_{m}(s) is obtainable. The general procedure of deriving T_{m}(s) is discussed in many textbooks, such as the one by BOYCE and DIPRIMA [8]. In what follows in this section, it is shown that a special finite element (i.e., TMB element) can be developed from Eq. (D3).
Let us consider a section of pulley, which is located in s_{0} ≤ s ≤ s_{L}, and let
T_{m}(S_{L})=  [  T_{UU}  T_{UF}  ]  (D5) 
T_{FU}  T_{FF}  
ψ_{m}(S_{L})=  [  U_{L}  ]  (D6)  
F_{int L}  
ψ_{m}(S_{0})=  [  U_{o}  ]  (D7)  
F_{int 0} 
P_{m} = T_{m}(s_{L}) 
s_{L} 
T_{m}(X)^{1 }B(x)dx (D8) 
Where
U_{L} = U_{m}(S_{L})  (D9) 
U_{o} = U_{m}(S_{0})  (D10) 
F_{int L} = F_{m}(S_{L})  (D11) 
F_{int 0} = F_{m}(S_{0})  (D12) 
Then Eq. (D3) can be rewritten as
[  U_{L}  ]  =  [  T_{UU}  T_{UF}  ][  F_{int L}  ]  + P_{m} (D13) 
F_{int L}  T_{FU}  T_{FF}  F_{int 0} 
Reorganizing Eq. (D13), we have
[  I  T_{UU}  ][  U_{L}  ]  =  [  0  T_{UF}  ][  F_{int L}  ]  + P_{m} (D14) 
0  T_{FU}  U_{0}  I  T_{FF}  F_{int 0} 
from which we can obtain
K_{m}D_{m} = F_{int m} + F_{ext m} (D15)
Where
K_{m} =  [  0  T_{UF}  ]^{1}[  I  T_{UU}  ]  (D16) 
I  T_{FF}  o  T_{FU}  
D_{m} = (U^{T}_{L}, U^{T}_{0})^{T} (D17)  
F_{int m} =  [  F_{int L}  ]  (D18)  
F_{int 0}  
F_{ext m} =  [  0  T_{UF}  ]^{1}  P_{m}  (D19)  
1  T_{FF} 
It is seen that Eq. (D15) is formulated in a finiteelement form, where K_{m} is element stiffness matrix, D_{m} is the element displacement vector, F_{int m} is the element internal force vector, and F_{ext m } is the element external force vector, The detailed procedures of developing TMB elements for shaft, enddisk and cylindrical shell are given in Appendices A, B and C, respectively.
Remarks:
 It is noted that in developing TMB finite elements the trial function for representing element displacement field is not required. This is the feature that distinguishes the MTM method from the conventional finite element method.
 In general, for any type of structures, whose governing differential equations can be reduced to the form of Eq. (D1), which can be further transformed into a set of constant ODEs, the corresponding transfer matrix T_{m}(s) of Eq. (D3) can be derived following some standard procedure described in [8]. In most cases, the procedure of deriving T_{m}(S) can be fulfilled by using a computer program, and P_{m} defined in Eq. (D8) can be evaluated analytically or numerically; and Therefore, a computer program can be made, where the input is the matrix H_{m} and the external loads, and the output is the TMB element stiffness matrix and element external force vector. In this program, there is no discretization error in computing TMB element stiffness matrix. As a FORTRAN subroutine, such a computer program has been developed in PSTRESS 3.0.