Transverse Crack Modeling and Validation in Rotor Systems , Including Thermal Effects

This article describes a model that allows the simulation of the static behavior of a transverse crack in a horizontal rotor under the action of weight and other possible static loads and the dynamic behavior of cracked rotating shaft. The crack breathes—that is, the mechanism of the crack’s opening and closing is ruled by the stress on the cracked section exerted by the external loads. In a rotor, the stresses are time-dependent and have a period equal to the period of rotation; thus, the crack periodically breathes. An original, simplified model allows cracks of various shapes to be modeled and thermal stresses to be taken into account, as they may influence the opening and closing mechanism. The proposed method was validated by using two criteria. First the crack’s breathing mechanism, simulated by the model, was compared with the results obtained by a nonlinear, threedimensional finite element model calculation, and a good agreement in the results was observed. Then the proposed model allowed the development of the equivalent cracked beam. The results of this model were compared with those obtained by the three-dimensional finite element model. Also in this case, there was a good agreement in the results. Therefore, the proposed models of the crack and the equivalent model of the beam can be inserted into the finite element model of the beam used for the rotor’s dynamicbehavior simulation; the obtained equations have timedependent coefficients, but they can be integrated into the frequency domain by using the harmonic balance method.

The description of the dynamical behavior of horizontal heavy rotors presenting a transverse crack has been presented by many researchers (see e.g., the overviews of Wauer (1990), Gasch (1993), and Dimarogonas (1996)), particularly as regards to the "breathing mechanism."This characteristic behavior takes place during the rotor rotation: the crack moves from the upper position in which the static bending moment forces the crack to be "closed" to the opposite position in which the crack is forced to be "open."The gradually opening and closing of the crack is called the breathing mechanism.Therefore, the stiffness changes periodically during one rotation, and correspondingly, so does the static deflection due to the weight (and to the bearing alignment conditions).In fact, the stiffness of the rotor with the open crack is smaller with respect to the closed crack situation, in which the stiffness equals the value of the uncracked shaft.The periodical deflection and stiffness can both be expanded in a Fourier series, of which the most important components are the 1x, 2x, and 3x rev.components.It can be shown that the same forces which excite the static deflection components also excite the dynamical vibrations when the shaft is rotating at higher speeds.
In addition to the effect of the static bending moments, thermal transients also have an influence on the breathing mechanism.Often, in real machines, the sudden change in vibrational behavior during thermal transients allowed, a crack in rotating shafts to be discovered.A case history is presented in Lapini et al. (1993) from which Figure 1 is derived.This shows the measured vibrations in a bearing of a cracked generator of a power plant, operating at rated speed, due to two thermal transients: first a

FIGURE 1
Cracked rotor thermal sensitivity.cooling in the refrigerating fluid from 40 • C to 34 • C, then a heating from 34 • C to 44 • C.During the heating transient (a positive temperature gradient), compressive stresses arise on the "skin" of the shaft and tensile stresses in the "core" of the shaft.This has the following effects: • If a local bow is presented in correspondence with the crack, due to the crack propagation, the thermal stresses reduce the existing bow and therefore a reduction in 1x rev.vibration component is obtained.
• The thermal stresses tend prevent and reduce the breathing behavior.Therefore a reduction of all components (1x, 2x, and 3x rev.) due to the breathing mechanism is generally obtained.
In the case of a cooling transient (a negative temperature gradient), the opposite situation with tensile stresses on the "skin" and compressive stresses on the "core" is obtained.This affects the crack in the following ways: • A local bow is generated (or the existing local bow is increased) and the 1x rev.component of the vibrations due to the unbalance and to the bow is modified.Previous considerations suggest that a reliable model for crack breathing mechanism that also takes into account the effect of thermal transients would prove useful.A possible model is described in this article.First, a simplified 1-D model, which was already presented in Bachschmid et al. (2000), is described to determine the temperature distribution in a cylindrical rotor.Then, the stiffness changes in the cracked element of the ro-

THERMAL BEHAVIOR AND BREATHING MECHANISM OF CRACKED ROTORS
To determine the temperature distribution, the equation of the thermal exchange is used in the case of axial-symmetry and of an infinite cylinder: where ρ is the density, c p the specific heat, k the conduction coefficient, r the radial coordinate, and T the generic temperature of the section.
To solve the equation, the finite differences have been used with two different boundary conditions depending on the type of heat exchange: a.In the case of convection on the external surface: where H is the thermal exchange coefficient between steam or surrounding fluid and metal, T rest is the temperature of the external surface of the section, and T w is the fluid temperature.b.In the case of an imposed temperature gradient T /T on the external surface: If the number of the discretization points in the finite difference method and the time step t are fixed and the material properties are considered as constant (at the average working temperature), the solving system is of the following type: where matrix [A] has constant coefficients and {T (t + t)} and {T (t)} are the vectors of the temperatures in the radial coordinates determined over the section.Note that the radial coordinates have been chosen in order to divide the section in to rings with equal surface.Therefore, they are thinner close to the external surface.By iterating Equation ( 4) with time step t, it is possible to determine the temperature distribution in each point of the section as a function of the time.The axial stress distribution corresponding to the temperature distribution is given by: where α is the linear thermal expansion, E the Young's modulus, ν the Poisson's coefficient, T (r ) the temperature at the radial coordinate r , and T m the average section temperature.The effect of the centrifugal forces on the axial stresses has been neglected.The thermal bow is calculated in the following way: The cracked areas in which no contact occurs cannot transmit stresses, therefore, the thermal stress distribution is no longer axial-symmetrical and the resulting moments M Bx and M By of this thermal stress distribution can be calculated with respect to the reference frame which is fixed to the beam.These moments, with opposite sign, are then applied to both ends of the equivalent beam of length l c , which will be defined later, and generate the bow.
This procedure is obviously a rough approximate approach, since the actual stress distribution over the cracked area is completely different from that assumed by the proposed simplified

FIGURE 2
Skin temperature (left) and axial stresses (right) for temperature gradient of ±100 • C/min.approach.However, it can be assumed that the bow depends on the uncompensated thermal stresses and the results confirm that the proposed procedure allows one to estimate the bow with satisfactory accuracy, as will be shown later.
The results obtained are different depending on short-time transients, in which the heating or cooling process involves more the outer skin of the beam than the whole body, or longtime transients, in which the complete body is interested by the temperature change and also high radial strains arise.With respect to the short transients, the stress distribution around the crack in long-lasting transients is different so as the deflections.

NUMERICAL SIMULATIONS
Some 1-D and 3-D calculations have been made on a round beam test specimen with a diameter of 25 mm, and a long presenting a crack of 50 mm with a relative depth of 25% and 50% of the diameter.Furthermore, mechanical loads have been applied to an extension (on the left-hand side) in order to avoid local deformations due to load application.The loads were a bending moment of 10 Nm and a torsion of 25 Nm.
Figure 2 shows the skin temperature and the corresponding axial stresses as a function of time for ±100 • C/min temperature gradient, calculated with the 1-D model, starting from a uniform temperature distribution of 40 • C.
Figure 3 shows the temperature and the stress distribution over the cross section of the beam at 5 s for the ±100 • C/min transient.
All  The behavior is as expected.During the negative transient, the positive tensile stresses on the skin and the negative compression stresses on the internal part vanish in correspondence of the crack.The crack is completely open.This can be seen in Figure 6 which shows the relative axial displacements on the crack surface: only in the dark blue area are the relative displacements vanishing small, so that some contact could occur.Figure 7 shows the axial stresses on the cracked section; on the cracked surface the axial stress is roughly 0, in the inner part a maximum compressive stresses of 11.8 MPa is reached, and on the outer part, close to the crack tip, a maximum of tensile stress of 117 MPa is reached.All these figures and some of the following have been obtained with the nonlinear 3-D model, for the crack of 50% depth.
During the positive transient, the negative compression stresses on the skin remain in correspondence with the crack, but, below that small portion, the rest of the crack has no stresses and is therefore open.Obviously, correspondence with the crack tip, very high tensile stress occur.
Figure 8 shows the stress distribution along the vertical diameter passing through the middle of the crack, in the case of positive and negative transients.This figure shows that the positive temperature transient effects the closure of the crack, as known from field experience.In fact, in a small area close to the skin, negative compression stresses appear.However, very high tensile stresses are generated at the crack tip (between 110 and 120 MPa), which could be responsible for the propagation of the crack in many machines (turbogenerators) in which a slowly propagating crack has been found.This then happens when the rotor is heated during the start-up procedure.
In the simplified model, 1-D thermal stresses are simply superposed to stresses arising from mechanical loads; the nonlinear effects cannot be accounted for.Therefore, excellent agreement is found in all parts which are not close to the crack (and to the boundary), but poor agreement is found on the cracked section.Despite this, acceptable agreement is found in the breathing mechanism, as shown hereafter, and a rather good agreement has been found in the evaluation of the bow due to smooth or sharp temperature transients (±20 • C/min and ±100 • C/min).
The deflections are calculated when the thermal transient is superposed to a mechanical load; this situation, which is always present in real machines, allows one to verify nonlinear effects in the superposition.In order to emphasize the thermal effect, the deflection due to the mechanical load alone has been subtracted.Only the position of external loads, which leads to open crack, has been considered.
From the results of Table 1, the following conclusions can be drawn.

FIGURE 6
Relative axial displacements on crack surface.

FIGURE 7
Axial stresses on cracked section.

FIGURE 8
Stress distribution along the vertical diameter passing though the middle of the crack.
• The 1-D linear model systematically underestimates the deflections for deeper cracks (50%) and for smaller cracks (25%) in the case of the sharper temperature gradients, • Since the errors introduced by the 1-D linear model are between a minimum of +0.7% and a maximum of 23%, the results can be considered a rough estimate of the actual behavior.
Other simulations show that when the cracked beam is without mechanical loads, loaded only by thermal stresses, then the symmetry of behavior is lost.Negative gradients produce the same deflections shown in Table 1, but positive gradients produce only very little deflections.
Thermal transients are considered long when the heat is propagated to the whole body of the beam.As said before, completely different stress and strain distributions can be found with respect to short transients.The axial stress distribution over the cross section of the beam, as obtained by the 1-D model with the smoother thermal transient (±20 • C/min) after 15 min, is very similar to the one obtained for ±100 • C/min gradient after 5 s and represented in Figure 3.
Considering the 50% crack and the smoother thermal heating transient (+20 • C/min) after 15 min, the contact pressures shown in Figure 9 are obtained, to be compared with Figure 10 related to the same transient, but after 10 s.Dark blue indicates no contact pressure.If the mechanical external load is superposed, then higher compressive stresses arise in the middle of the beam (the distribution over the cracked surface is shown in Figure 11).The behavior is therefore as would be expected.

THE 1-D CRACK MODEL
The proposed crack model is briefly described.The different steps for modeling the breathing behavior, including thermal effects, are the following: 1.The bending moment M, due to the weight and the bearing alignment conditions of the rotor, is calculated in correspondence to the cracked section.angular position ϑ of the main axis of inertia can be found.2.5 Now the procedure from 2.1 to 2.5 is repeated with the new value of ϑ, until ϑ converges to a stable value.3.At this point, the position of the main axis, the second moments of area and M Bx , M By , are known.The second moment of area J x , J y , and J xy with respect to the fixed reference frame (xy) and the components of the moments due to the thermal stress distribution (M x , M y ) with respect to the same reference frame, are calculated.This will be repeated for each angular position of the shaft.
The second moments of area, which are a function of the angular position, can then be used for calculating the stiffness K c ( t) that is also a function of the angular position of the cracked beam element, which has a suitable length l c .The moments of the thermal stress distribution are used for calculating the thermal bow.

EQUIVALENT BEAM STIFFNESS MATRIX
Once the breathing mechanism and the second moments of area have been defined for the different angular positions, the stiffness matrix of the cracked element of suitable length l c can be calculated, assuming a Timoshenko beam.

FIGURE 11
Contact pressure distribution over the cracked section, transient +20 • C/min, 10 s.

FIGURE 12
Cracked cross section.
The stiffness matrix (square, symmetrical, 12 × 12 elements) is represented in Equation ( 6): where the coefficients are defined as: The parameter φ accounts for the shear effects and is given by: E and G, respectively, are the Young's modulus and the shear modulus and S is the cross section area.The different lengths l c , l a , and l t responsible for the direct stiffnesses and the parameters c 1 , c 2 , and c 3 responsible for the cross-coupling terms, have been tuned by means of the 3-D model.The lengths l c , l a , and l t in some kind represent the extent to which the crack exerts his influence.The lengths are obviously related to the depth of the crack as shown, e.g., in Figure 13, where the relative length (ratio of length l c to diameter) is represented as a function of the relative depth (ratio of depth to diameter).The cracked element can then be introduced in the original beam element of length l.If the length l is greater than the maximum value among l c , l a , and l t , the stiffness of the element is easily calculated as also its deflection due to external loads in the different angular positions.If, on the contrary, the original element is shorter than the same maximum value, some additional processing of K c ( t) is required, as shown below.

DERIVATION OF LOCAL CRACK STIFFNESS
The rotation-dependent stiffness matrix K c of a cracked beam element of length l (Figure 14a) can be split into a stiffness matrix with three different parts: the stiffness of two equal uncracked beams of length l/2 combined with a "local" stiffness part composed by "springs" which represent the crack (Figure 14b).

FIGURE 13
Equivalent length l c of the cracked beam.
These springs have an infinite stiffness when the crack is closed and a finite stiffness when the crack is open.If we want to extract the local stiffness due to the crack from total stiffness, we have to introduce two additional nodes, namely nodes L and R, representing, respectively, the left-and right-hand sides of the crack faces.
Denoting by 1 the initial node of the cracked beam element, by 2 its final node, by X 1 and X 2 the vectors of the displacements of node 1 and 2, respectively, we can consider K c a 12×12 matrix which leads to following equation: where F 1 and F 2 represent the force vectors applied to the end nodes of the beam.This matrix, being 1 and 2 external nodes of the beam, has to be distinguished from local crack stiffness K II (see Figure 14b) corresponding to nodes L and R: Considering the stiffness matrix according to (Figure 14b) we get the following structure: This matrix can be reordered considering the external (E) and internal (I) degrees of freedom.Calling we get: where K EI and K IE are the stiffness matrices of the uncracked beam elements of length l/2.Therefore K II can be extracted according to the following procedure: combining Equation ( 13) with Equation ( 17) from which we extract K II : It is interesting to note that the resulting local stiffness matrix K II does not show coupling terms between rotations and deflections, which instead exist in matrix K c .This result was expected, considering that the local stiffnesses refer to a "beam" with zero length.
The local stiffness matrix can then be introduced in the cracked beam element remaining any restriction about the length of the element.The bow which has been calculated with the equivalent length l c beam generates a relative angular deflection

BREATHING MECHANISM MODEL VALIDATION WITH TEMPERATURE GRADIENTS
Positive and negative gradients have been considered, and their effects on 25 and 50% depth cracks.Since the dimensions of the model are small, and the maximum statical bending stress

FIGURE 15d
50% depth crack, after 5 s, temperature transient on the outer surface −100 • C/min, with external loads, angle 180 • .reach 6.5 N/mm 2 , a rather robust thermal transient has been applied, in order to generate stresses of the same order of magnitude.
Constant gradients of ±100 • C/min and ±20 • C/min have been applied to the external surface of the cylinder and the resulting stresses and the breathing behavior have been calculated after 5 and 10 s, respectively.
If the external load (bending + torsion) is superposed to the thermal load, the following figure is obtained for the 50% depth crack and the sharper temperature transient (−100 • C/min after 5 s): Figure 15 shows the situation corresponding to the following angular positions of 90 • , 120 • , 150 • , and 180 • .
For smaller angles, the crack is completely open; for angles from 180 to 360 • a symmetrical behavior is found, as expected.Similar results are obtained for other crack depths and other thermal transients.Note that the dark blue zones indicate zones in which contact between the crack faces occurs.
Although the agreement is not so good as in the case of only mechanical loading, it can be considered satisfactory, bearing in mind that the actual stress distribution over the cracked section deviates strongly from superposition of stresses used by the simplified linear model.

CONCLUSIONS
The thermal behavior of cracked round beams has been analyzed by means of a 3-D nonlinear model.The temperature distribution is unaffected by the crack, but the stress and strain distributions are strongly influenced by the crack.This results in a bow of the beam.
A 1-D model for calculating temperature and axial stress distributions in an infinite cylinder is presented and the results of the axial stress distribution are used for calculating thermal bows, by means of a rough approximation.The acceptable agreement with results obtained by means of the 3-D nonlinear model, validates the proposed procedure.
Furthermore, a 1-D crack model has been described that allows cracks of different shapes to be modeled and thermal stresses to be taken into account.The breathing mechanism obtained with this 1-D model has been validated with 3-D nonlinear calculations.The obtained good agreement allows one to propose the 1-D model as a powerful tool for cracked rotor analysis.

ACKNOWLEDGMENTS
This work is partially funded by the MURST (Italian Ministry for the University and Scientific Research) Cofinanziamento "Identificazione di Malfunzionamenti in Sistemi Meccanici" for the year 1999.

NOMENCLATURE
tor are described by introducing a 6 d.o.f.(degrees of freedom) per node model.The proposed models are validated by means of the comparison of the results obtained with 3-D f.e.(finite elemeents) calculations.
these results have been obtained with the 1-D model, but the actual temperature and stress distributions, obtained by the 3-D model in cross sections that are unaffected by the crack and by the boundary conditions are quite similar.Figures 4 and 5 represent the axial stress distribution on the longitudinal section of the cracked beam, clamped on the righthand side end: Figure 4 in the case of negative sharp transient (−100 • C/min in 5 s) and Figure 5 in the case of positive sharp transient (+100 • C/min in 5 s).

FIGURE 3
FIGURE 3Internal temperatures (left) and axial stresses (right) for ±100 • C/min gradient at time 5 s.

•
The thermal stresses tend to hold the crack open.The breathing mechanism is modified and certainly an increase in 2x rev.component is expected.

TABLE 1
Comparison of deflections due to thermal stresses in[µm] c 1 , c 2 , c 3 cross-coupling stiffness coefficient tuning para-length of the equivalent beam used to model the crack l t equivalent length for the torsional stiffness M bending moments due to the weight and bearing alignment conditions of the rotor M x , M y bending moments due to thermal stress distribution, with respect to the fixed reference M Bx , M By bending moments due to thermal stress distribution, with respect to the main axes , Y g coordinates of the center of area of the section X gT , Y gT coordinates of the center of area of the section that participates to the torsion