# Finite Element Multibody Simulation of a Breathing Crack in a Rotor with a Cohesive Zone Model

Authors: Rugerri Toni Liong; Carsten Proppe
Journal: ISRN Mechanical Engineering
 (2013)

---

## Abstract

The breathing mechanism of a transversely cracked shaft and its influence on a rotor system that appears due to shaft weight and inertia forces is studied. The presence of a crack reduces the stiffness of the rotor system and introduces a stiffness variation during the revolution of the shaft. Here, 3D finite element (FE) model and multibody simulation (MBS) are introduced to predict and to analyse the breathing mechanism on a transverse cracked shaft. It is based on a cohesive zone model (CZM) instead of linear-elastic fracture mechanics (LEFM). First, the elastic cracked shaft is modelled by 3D FE. As a second step, the 3D FE model of the shaft is transferred into an MBS model in order to analyze the dynamic loads, due to the crack, and the inertia force acting during rotation at different rotating speeds. Finally, the vibration responses in the centroid of the shaft obtained from MBS have been exported into FE model in order to observe the breathing mechanism. A bilinear crack closure model is proposed. The accuracy of the bilinear crack closure model and the solution techniques have been demonstrated by a comparison with the corresponding results of previous publications.


---

## Body

## 1. Introduction

Fatigue cracking of rotor shafts has long been identified as a limiting factor for safe and reliable operation of turbomachines. It can lead to catastrophic failure and great economic loss if not detected early. A crack in the rotor causes local changes in stiffness. These changes, in turn, affect the dynamics of the system: frequency of the natural vibrations and the amplitudes of forced vibrations are changed. If a cracked shaft rotates under external loading, the crack opens and closes regularly during the revolution of the shaft; it breathes. The breathing mechanism is produced by the stress distribution around the crack mainly due to the action of bending moment, while the effect of torsion is negligible. Usually, shaft cracks breathe when crack sizes are small, running speeds are low, and radial forces are large [1].The influence of a breathing crack on the vibration of a rotating shaft has been in the focus of many researchers. Comprehensive literature survey of various crack modelling techniques, system behaviour of cracked rotor, and detection procedures to diagnose fracture damage were contributed Wauer [2], Dimarogonas [3], Sabnavis et al. [4], and Kumar and Rastogi [5]. More recent studies have been reviewed by Bachschmid et al. [6]. They noted that the breathing mechanism of cracks in rotating shafts can be accurately investigated by means of 3D nonlinear finite element models. Based on these simulations, simple approximation of the breathing mechanism, describing the location and extension of the crack closure line during rotation, can be established.The breathing mechanism is the result of the stress and strain distribution around the cracked area, which is due to static loads, like the weight, bearing reaction forces, and dynamical loads, and due to the unbalance- and the vibration-induced inertia force distribution [7]. The accurate modelling of the breathing behaviour, that is, the gradually opening and closing of the crack during one rotation of the shaft, still needs further investigations. An original method for calculating the constitutive law of a cracked beam section under bending has been proposed by Andrieux and Varé [8]. Based on three-dimensional computations taking into account the unilateral contact between the lips of the crack, it consists in defining a constitutive relation between the bending moment applied to the cracked section and the resulting field of displacements, compatible with the beam theory so that it can be used in rotor-dynamics software. Varé and Andrieux [9] extended this method in order to show how shear effects can be implemented in the model of a cracked section. Arem and Maitournam [10] presented stiffness variations deduced from 3D FE calculations accounting for the unilateral contact between the crack lips as originally developed by Andrieux and Varé [8]. Dimarogonas and Paipetis [11] introduced the FE models of the rotating shaft cracks. They used the fracture mechanics and obtained a full 6×6 flexibility matrix for a transverse open surface crack on shaft. Darpe et al. [12] studied the coupling between longitudinal, lateral, and torsional vibrations together for a rotating cracked shaft. By including the axial degree of freedom in their analysis, the stiffness matrix formulated is an extension of the one developed by Ostachowicz and Krawczuk [13]. Kulesza and Sawicki [14] used the rigid FEM for modelling a crack in a rotating shaft. The crack is presented as a set of spring-damping elements of variable stiffness connecting two sections of the shaft. Bouboulas and Anifantis [15] develop a finite element model in order to study the vibrational behavior of a beam with a nonpropagating edge crack. According to their model, the beam is discretized into finite elements while the breathing crack behavior is treated as a full frictional contact problem between the crack surfaces.The cohesive zone model (CZM) has been widely used as an alternative to stress intensity factor-based fracture mechanics. It can deal with the nonlinear zone ahead of the crack tip due to plasticity or microcracking. The CZM describes material failure on a more phenomenological basis (i.e., without considering the material microstructure). The general advantage of this model when compared to classical fracture mechanics is that the parameters of the respective models depend only on the material and not on the geometry. This concept guarantees transferability from specimen to structure over a wide range of geometries. The origin of the cohesive zone concept can be traced back to the strip yield model proposed by Dugdale [16] and Barenblatt [17] in which the narrow zone of localized deformation ahead of the crack tip was substituted by cohesive traction between the bounding surfaces. The constitutive behaviour which causes the cohesive elements to open and eventually to fail is described by the so-called traction separation law. It relates the traction vector to the displacement jump across the interface and is usually called separation. The energy dissipated by the element until total failure is derived as the integral of the traction-separation. However, the traction-separation law depends on the stress state, which can be characterised by the triaxiality (i.e., the hydrostatic stress divided by the von Mises equivalent stress). This issue was first investigated by Siegmund and Brocks [18, 19]. The approach was extended to simulation of dynamic ductile crack growth by Anvari et al. [20], Scheider [21]. Banerjee and Manivasagam [22] proposed a versatile CZM to predict ductile fracture at different states of stress. The formulation developed for mode-I plane strain accounts explicitly for triaxiality of the stress state by using basic elastic-plastic constitutive relations combined with two new model parameters, which are independent of the stress state.The theory of strain energy release rate and SIF combined with rotordynamics built the foundation of the dynamic of cracked rotors based on LEFM [23]. This theory has two major limitations; namely, it can only be used if there is an initial crack, and the fracture process zone must be small compared with the dimensions of the shaft. Due to geometrical complexity, some simplifications had been made for the crack profile, such as straight-edged, circular and elliptical crack model to analyze such crack problem. Since analytical SIFs for an edge crack in a rotating cylinder are not available, the shaft is considered to be a sum of elementary independent rectangular strips (Andrieux and Varé  [8]) and no interaction between them is assumed to take place (Chasalevris and Papadopoulos [24]). The geometric functions that describe the strain energy density are often not accurate enough, due to the fact that the crack passes from stress state caused by the vertical moment to that of a horizontal moment. Then, the compliance is obtained by integrating along the crack tip. If the crack depth exceeds the radius of the shaft, then the elements of the compliance matrix present a divergence. This is due to the singularity that the strain energy release rate method has near the edges of the crack tip, giving thus the infinite values. It was reported by Papadopoulos [25] that divergence does not reflect reality. The crack tip is supposed to be formed by the boundary between the cracked areas and the uncracked areas for the regions in which the breathing crack is open. However, the SIF will not appear at the boundary between the closed cracked areas and the open cracked areas [7]. Liong and Proppe [26] proposed a method for the evaluation of the stiffness losses in a rotor with a transverse breathing crack. Their method is based on a CZM and accounts explicitly for triaxiality of the stress state by using constitutive relations.The aim of the present study is to propose a method for the evaluation of stiffness losses of the cracked shaft based on a CZM. The method provides a direct relationship between the material properties of the shaft and the relative crack depth, thus representing the physical phenomenon closely. This paper introduces a method based on a 3D FE model and MBS in order to study the breathing mechanism of an elastic cracked shaft. The CZM formulation is implemented in a 3D FE model. The elastic cracked shaft with various relative crack depths is modelled by 3D FE and then transferred into an MBS model in order to analyze the dynamic loads acting during rotation at different rotating speeds. The vibration responses in the centroid of the rotating shaft obtained from MBS have been exported into FE to observe the breathing mechanism. The accuracy of the results is demonstrated through comparisons with the results available in the literature.


## 2. The Cohesive Zone Model

### 2.1. Basis of the Cohesive Process Zone Model

In the CZM, fracture nucleates as discontinuity surface able to transmit tensile load before opening above a given displacement. Formation and extension of this surface require that the maximum principal stress reaches a given value, namely, the cohesive strength of the material. When this occurs, the crack surface initiates or grows perpendicularly to the direction of the maximum principal stress. The two faces of the surface exert on each other equal and opposite tensile stresses (cohesive stresses) whose value is a unique functionf(δ) of the separation δ between the faces. When the separation reaches another given value (the critical separation, δc), the cohesive stress vanishes and fracture takes place. Fracture consists of the initiation and propagation of a crack produced by the opening and advance of the cohesive zone (the zone where the cohesive stresses act) ahead of the crack tip as shown in Figure 1.Figure 1
Cohesive crack process zone.The fracture behaviour of each material is described by the cohesive traction as a function(1)σn=σmaxf(δ),
where σmax is the peak value of traction. The function f(δ) defines the shape of the traction-separation law and the area under the cohesive law curve is the work of separation or cohesive energy Gc:
(2)Gc=∫0δc‍σn(δ)dδ.Ideally, the model within the cohesive zone should be able to replicate the constitutive behaviour of the undamaged material: linear elastic followed by strain hardening, till the conditions for the initiation of softening due to damage are reached. The softening process representative of increasing material degradation is triggered by rapid growth of voids as consequence of a highly triaxial state of stress.Since the CZM is a phenomenological model, various formulations for defining the shape of traction-separation law and the cohesive values are in use [27]. A versatile CZM to predict ductile fracture at different states of stress is proposed in [22]. The formulation developed for mode-I plane strain accounts explicitly for triaxiality of the stress state by using basic elastic-plastic constitutive relations combined with two stress-state-independent new model parameters. The proposed traction-separation law has three distinct regions of constitutive behaviour: the traction separation law is linear up to the separation limit δn=δ1 and exhibits strain hardening up to δn=δ2 followed by a softening curve. The relevant variables and zones are sketched in Figure 2:
(3)σn=(1+3χef)2E3δfor0<δ<δ1.Figure 2
Traction separation law for ductile materials.The linear behaviour of the cohesive zone exists till the separation limit defined by the von Mises yield condition is reached:(4)δ1=3σY2E.Further separation results in strain hardening up to(5)δ2=32(Ce-(3/2)χef+σYE).


### 2.2. Numerical Representation of Cohesive Zone Models

When fracture proceeds, energy must be supplied by external loads. The bounding material undergoes elastoplastic deformation involving elastic energy and plastic dissipative energy. In addition to plasticity, energy is supplied to the fracture process zone in the form of cohesive energy that is dissipated within the cohesive elements. The cohesive energy is the sum of the surface energy and all dissipative processes that take place within the crack tip regime. For the present problem, a perfect energy balance between external workW and the sum of elastic energy Eel and cohesive energy Gc [27, 28] will be assumed. The energy balance is given by
(6)W=Eel+Gc,(7)W=∫0t(∫Vσ:ε˙eldV)dt+∫0t(∫Sσn:δ˙dS)dt,
where σ, ε˙el, σn, and δ˙ are nominal stress tensor, elastic strain rate, cohesive traction, and cohesive separation rate, respectively (including the terms for specimen volume V and the internal specimen surface S). The external work due to the applied force is given by
(8)W=∫0t(∫St·vdS)dt+∫0t(∫Vb·vdV)dt,
where  v  is the velocity field vector, t the exterior surface traction vector, and b the body force vector.The cohesive surface contribution, representing the crack and the process zone in front of the crack tip, is described by the integral over the internal surfaceS. In this formulation, σn denotes the cohesive strength, that is, the maximum traction value that can be sustained within the cohesive zone. The cohesive length δ is the value of the displacement jump across the crack surfaces at which the stress carrying capacity of the cohesive elements reaches its maximum value. By creating new surfaces, the traction and the stiffness of the cohesive zone elements connecting these newly created surfaces are made to vanish, but the displacements across them are still continuous. During FE analysis, the amount of external work, elastic work, and cohesive energy is calculated. The energy balance given by (7) and (8) is maintained in all FE computations.


### 2.3. Finite Element Implementation

The cohesive surface contribution as shown in Figure3 is implemented into an FE code for eight node elements based on CZM concepts. The fine mesh around the crack tip is shown in Figure 4. After making convergence studies, the shaft is discretized with 4108 plane strains and 8 node quadrilateral elements. The length of elements is ℓ, where the length of shaft is 20ℓ. Near the crack tip, the size of the element is ℓ/20. A total of 4833 nodes are used to model the geometry as shown in Figure 4. The fracture process zone is modeled by 8 node rectangular cohesive elements having zero thickness. 38 cohesive elements are used. One face of cohesive elements is connected to continuum elements with 4 nodes, while the other face is given symmetric displacement boundary conditions. Thus, an artificial interface is created along the fracture process zone.Figure 3
Representation of the fracture process using CZM in FE.Figure 4
Finite element mesh model near the crack tip.The maximum cohesive traction is assumed to be the same as the yield strength of the material of shaft; namely,σn=250 MPa. The maximum separation at the end of the elastic zone is  δ1=10 μm. The interface stiffness or penalty stiffness Kp can be obtained directly from  Kp=σn/δ1  and is used to ensure a stiff connection between the surfaces of the material discontinuity. The values of the penalty stiffness should be selected in a way such that the connection between the two elements and numerical stability are guaranteed.


## 2.1. Basis of the Cohesive Process Zone Model

In the CZM, fracture nucleates as discontinuity surface able to transmit tensile load before opening above a given displacement. Formation and extension of this surface require that the maximum principal stress reaches a given value, namely, the cohesive strength of the material. When this occurs, the crack surface initiates or grows perpendicularly to the direction of the maximum principal stress. The two faces of the surface exert on each other equal and opposite tensile stresses (cohesive stresses) whose value is a unique functionf(δ) of the separation δ between the faces. When the separation reaches another given value (the critical separation, δc), the cohesive stress vanishes and fracture takes place. Fracture consists of the initiation and propagation of a crack produced by the opening and advance of the cohesive zone (the zone where the cohesive stresses act) ahead of the crack tip as shown in Figure 1.Figure 1
Cohesive crack process zone.The fracture behaviour of each material is described by the cohesive traction as a function(1)σn=σmaxf(δ),
where σmax is the peak value of traction. The function f(δ) defines the shape of the traction-separation law and the area under the cohesive law curve is the work of separation or cohesive energy Gc:
(2)Gc=∫0δc‍σn(δ)dδ.Ideally, the model within the cohesive zone should be able to replicate the constitutive behaviour of the undamaged material: linear elastic followed by strain hardening, till the conditions for the initiation of softening due to damage are reached. The softening process representative of increasing material degradation is triggered by rapid growth of voids as consequence of a highly triaxial state of stress.Since the CZM is a phenomenological model, various formulations for defining the shape of traction-separation law and the cohesive values are in use [27]. A versatile CZM to predict ductile fracture at different states of stress is proposed in [22]. The formulation developed for mode-I plane strain accounts explicitly for triaxiality of the stress state by using basic elastic-plastic constitutive relations combined with two stress-state-independent new model parameters. The proposed traction-separation law has three distinct regions of constitutive behaviour: the traction separation law is linear up to the separation limit δn=δ1 and exhibits strain hardening up to δn=δ2 followed by a softening curve. The relevant variables and zones are sketched in Figure 2:
(3)σn=(1+3χef)2E3δfor0<δ<δ1.Figure 2
Traction separation law for ductile materials.The linear behaviour of the cohesive zone exists till the separation limit defined by the von Mises yield condition is reached:(4)δ1=3σY2E.Further separation results in strain hardening up to(5)δ2=32(Ce-(3/2)χef+σYE).


## 2.2. Numerical Representation of Cohesive Zone Models

When fracture proceeds, energy must be supplied by external loads. The bounding material undergoes elastoplastic deformation involving elastic energy and plastic dissipative energy. In addition to plasticity, energy is supplied to the fracture process zone in the form of cohesive energy that is dissipated within the cohesive elements. The cohesive energy is the sum of the surface energy and all dissipative processes that take place within the crack tip regime. For the present problem, a perfect energy balance between external workW and the sum of elastic energy Eel and cohesive energy Gc [27, 28] will be assumed. The energy balance is given by
(6)W=Eel+Gc,(7)W=∫0t(∫Vσ:ε˙eldV)dt+∫0t(∫Sσn:δ˙dS)dt,
where σ, ε˙el, σn, and δ˙ are nominal stress tensor, elastic strain rate, cohesive traction, and cohesive separation rate, respectively (including the terms for specimen volume V and the internal specimen surface S). The external work due to the applied force is given by
(8)W=∫0t(∫St·vdS)dt+∫0t(∫Vb·vdV)dt,
where  v  is the velocity field vector, t the exterior surface traction vector, and b the body force vector.The cohesive surface contribution, representing the crack and the process zone in front of the crack tip, is described by the integral over the internal surfaceS. In this formulation, σn denotes the cohesive strength, that is, the maximum traction value that can be sustained within the cohesive zone. The cohesive length δ is the value of the displacement jump across the crack surfaces at which the stress carrying capacity of the cohesive elements reaches its maximum value. By creating new surfaces, the traction and the stiffness of the cohesive zone elements connecting these newly created surfaces are made to vanish, but the displacements across them are still continuous. During FE analysis, the amount of external work, elastic work, and cohesive energy is calculated. The energy balance given by (7) and (8) is maintained in all FE computations.


## 2.3. Finite Element Implementation

The cohesive surface contribution as shown in Figure3 is implemented into an FE code for eight node elements based on CZM concepts. The fine mesh around the crack tip is shown in Figure 4. After making convergence studies, the shaft is discretized with 4108 plane strains and 8 node quadrilateral elements. The length of elements is ℓ, where the length of shaft is 20ℓ. Near the crack tip, the size of the element is ℓ/20. A total of 4833 nodes are used to model the geometry as shown in Figure 4. The fracture process zone is modeled by 8 node rectangular cohesive elements having zero thickness. 38 cohesive elements are used. One face of cohesive elements is connected to continuum elements with 4 nodes, while the other face is given symmetric displacement boundary conditions. Thus, an artificial interface is created along the fracture process zone.Figure 3
Representation of the fracture process using CZM in FE.Figure 4
Finite element mesh model near the crack tip.The maximum cohesive traction is assumed to be the same as the yield strength of the material of shaft; namely,σn=250 MPa. The maximum separation at the end of the elastic zone is  δ1=10 μm. The interface stiffness or penalty stiffness Kp can be obtained directly from  Kp=σn/δ1  and is used to ensure a stiff connection between the surfaces of the material discontinuity. The values of the penalty stiffness should be selected in a way such that the connection between the two elements and numerical stability are guaranteed.


## 3. Investigation of Breathing Mechanism

### 3.1. Method and Assumptions

The breathing mechanism of a simple cracked rotor on rigid supports has been investigated. The method of the investigation is summarized in Figure5. Several assumptions are made: the location of the crack is assumed to be known at the mid-span of the shaft which is the same as the maximum deflection point. Only a single transverse crack has been considered in the shaft. The plane strain condition is assumed at the crack front due to the geometry constraint. Both ends of the shaft are rigidly supported. The shaft has a uniform circular cross section. Length ℓ and diameter d of the shaft are assumed 1.0 m and 0.08 m, respectively. The material of the shaft is considered to be homogeneous and isotropic. Young’s modulus E, Poisson’s ratio ν, and mass density ρ are 210 GPa, 0.3, and 7850 kg/m3, respectively. Yield strength σY and ultimate strength of material σmax are 250 MPa and 400 MPa, respectively. The cohesive elements along the crack surfaces are implemented as interface elements as shown in Figure 6.Figure 5
Principle of the breathing crack simulation using integrated FE and MBS.Figure 6
FE model of an elastic cracked shaft.


### 3.2. Breathing Mechanism Simulation

The breathing mechanism generated by the rotating bending load used in the literature has some limitations. Due to the presence of inertia forces, the dynamic behaviour of rotating structures is different from those of static structures. Although in case of weight dominance, the amplitude of the vibration response due to inertia forces is smaller than the amplitude due to weight forces, elastic forces, and presence of a crack (if the shaft is assumed to be balanced), using the inertia force taken into account will yield more accurate results. The 3D FE calculations allow the breathing mechanism to be predicted accurately. The breathing mechanism is strongly influenced by the weight of the shaft [29]. In many publications on the breathing crack, assumptions on the breathing mechanism have been made. For transverse cracks in rotating shafts, the variation of the crack front line perpendicular to the crack front has been extensively addressed by Darpe et al. [12]. This breathing crack form is also used by Jun et al. [30] and Sinou and Lees [31]. An elliptical shape has been proposed by Shih and Chen [32]. Based on FE model simulation and some reported experimental results, the breathing crack shape is modelled by a parabolic shape that opens and closes due to bending stresses (Liong and Proppe [26, 33, 34]).The idea is that the vibration responses in the centroid of the shaft obtained from MBS are exported into FE model in order to analyse the breathing mechanism, as schematically shown in Figure7. The opening crack is simulated for one cycle of revolution of the cracked shaft specimen in steady-state condition. The breathing mechanism is generated by the bending due to external load (weight) by increasing the angle by steps of π/12 rad. The breathing (open and closed crack areas are evaluated in each angular step) is observed by the nodal displacement and the stress distribution (tensile or compressive stress) around the crack. The prediction of the breathing mechanism was performed by the following steps: for each angle of rotation, stress distribution due to the bending moment is recorded over the cross section. Compressive (negative) stresses indicate closed region. The crack opens where zero or very small positive numerical values of stresses appear and the contact forces vanish. Displacements and stresses can be observed at the crack surfaces. In order to avoid local deformations due to the application of loads, the deflection and stress distribution at each crack element area must be evaluated as shown in Figure 8.Figure 7
The transfer of the vibration responses from MBS into FE.Figure 8
Evaluation of displacements and stresses for state of crack.


## 3.1. Method and Assumptions

The breathing mechanism of a simple cracked rotor on rigid supports has been investigated. The method of the investigation is summarized in Figure5. Several assumptions are made: the location of the crack is assumed to be known at the mid-span of the shaft which is the same as the maximum deflection point. Only a single transverse crack has been considered in the shaft. The plane strain condition is assumed at the crack front due to the geometry constraint. Both ends of the shaft are rigidly supported. The shaft has a uniform circular cross section. Length ℓ and diameter d of the shaft are assumed 1.0 m and 0.08 m, respectively. The material of the shaft is considered to be homogeneous and isotropic. Young’s modulus E, Poisson’s ratio ν, and mass density ρ are 210 GPa, 0.3, and 7850 kg/m3, respectively. Yield strength σY and ultimate strength of material σmax are 250 MPa and 400 MPa, respectively. The cohesive elements along the crack surfaces are implemented as interface elements as shown in Figure 6.Figure 5
Principle of the breathing crack simulation using integrated FE and MBS.Figure 6
FE model of an elastic cracked shaft.


## 3.2. Breathing Mechanism Simulation

The breathing mechanism generated by the rotating bending load used in the literature has some limitations. Due to the presence of inertia forces, the dynamic behaviour of rotating structures is different from those of static structures. Although in case of weight dominance, the amplitude of the vibration response due to inertia forces is smaller than the amplitude due to weight forces, elastic forces, and presence of a crack (if the shaft is assumed to be balanced), using the inertia force taken into account will yield more accurate results. The 3D FE calculations allow the breathing mechanism to be predicted accurately. The breathing mechanism is strongly influenced by the weight of the shaft [29]. In many publications on the breathing crack, assumptions on the breathing mechanism have been made. For transverse cracks in rotating shafts, the variation of the crack front line perpendicular to the crack front has been extensively addressed by Darpe et al. [12]. This breathing crack form is also used by Jun et al. [30] and Sinou and Lees [31]. An elliptical shape has been proposed by Shih and Chen [32]. Based on FE model simulation and some reported experimental results, the breathing crack shape is modelled by a parabolic shape that opens and closes due to bending stresses (Liong and Proppe [26, 33, 34]).The idea is that the vibration responses in the centroid of the shaft obtained from MBS are exported into FE model in order to analyse the breathing mechanism, as schematically shown in Figure7. The opening crack is simulated for one cycle of revolution of the cracked shaft specimen in steady-state condition. The breathing mechanism is generated by the bending due to external load (weight) by increasing the angle by steps of π/12 rad. The breathing (open and closed crack areas are evaluated in each angular step) is observed by the nodal displacement and the stress distribution (tensile or compressive stress) around the crack. The prediction of the breathing mechanism was performed by the following steps: for each angle of rotation, stress distribution due to the bending moment is recorded over the cross section. Compressive (negative) stresses indicate closed region. The crack opens where zero or very small positive numerical values of stresses appear and the contact forces vanish. Displacements and stresses can be observed at the crack surfaces. In order to avoid local deformations due to the application of loads, the deflection and stress distribution at each crack element area must be evaluated as shown in Figure 8.Figure 7
The transfer of the vibration responses from MBS into FE.Figure 8
Evaluation of displacements and stresses for state of crack.


## 4. Results and Discussion

### 4.1. Breathing Mechanism Results

Figure9 displays the result of breathing mechanism during half-revolution of a rotating shaft (from closed to open crack) with relative crack depth  a/d=0.1. As can be seen the crack will start to appear and begin to open at π/6. The crack opens more slowly at the beginning, but increases its opening at π/3. At 5π/6 it is already completely open. Similarly, the crack closes again at 7π/6 and increases its closing at 3π/2. The crack is already completely closed at 11π/6. It is interesting to note that the angular positions where the crack opens again (from fully closed) are the same as where the crack begins to close.Figure 9
Breathing crack of the rotating shaft: relative crack deptha/d=0.1 at 600 rpm.Results of the relative crack depth versus the shaft rotation angle during breathing are shown in Figure10. Four functions are used to approximate the relative crack depth, namely, linear, quadratic, harmonic, and power functions. Table 1 gives the approximations for a relative crack depth a/d as a function of shaft rotation angle θ. It is shown that harmonic and quadratic functions have the minimal error and the nearest correlation with the simulation results. In this sense, the relative crack depth during crack opening should be understood not as linear increasing but as harmonic function or quadratic polynomial.Table 1
Curve fitting of the simulation results fora/d=0.1, θ in [rad].
Model
Mathematical model
SSE
R-square
Linear function
a
/
d = 0.039θ + 0.001
0.000 580
0.951 6
Quadratic polynomial
a
/
d = −0.012θ2 + 0.076θ − 0.020
0.000 011
0.999 0
Sinus function
a
/
d = 0.010sin (0.593θ − 0.145)
0.000 006
0.999 5
Power function
a
/
d = 0.123θ0.373 − 0.077
0.000 099
0.991 7Curve fitting of breathing mechanism results.
(a)
(b)
(c)
(d)


### 4.2. Discussion

Due to the presence of gravity, the upper portion of the cracked rotor at the beginning of a rotation is under compression and the crack is closed. As the rotor continues to rotate and the gravity direction is constant, the upper part now comes in the lower tensile region causing the crack to open. Moreover, the accuracy of the results is modelled by a bilinear model as shown in Figure11. The crack begins opening at shaft rotation angle θ0 and at crack front angle α that defines the angle between two crack front lines. The shaft rotates with shaft rotation angle θ and the crack front angle α increases until maximum crack width bmax is reached. The crack front angle α and the crack width b are nonlinearly increased during rotation from θ0 to θmax (Figures 12 and 13). One can obtain the crack front angle α and the crack width b as functions of the shaft rotation angle θ, respectively, by
(9)α={0,θ≤θ0crackstillclosed,-0.018θ3+0.016θ2+0.35θ+2.4,θ0≤θ≤θmaxcrackopens,0,θ≥θmaxcrack opens completely,bbmax={0,θ≤θ0crack still closed,0.076θ4-0.54θ3+1.1θ2-0.23θ-0.0056,θ0≤θ≤θmaxcrack opens,0,θ≥θmaxcrack opens completely.Figure 11
Bilinear crack closure model.Figure 12
Crack front angleα during rotation of shaft.Figure 13
Normalized crack widthb/bmax during rotation of shaft.For the relative crack deptha/d=0.1, the crack begins to open at  θ0=π/6  and  α=5π/6, while the crack opens completely at  θmax=5π/6  and  α=π. Therefore, the crack front angle and the crack width during crack opening are
(10)α=-0.018θ3+0.016θ2+0.35θ+2.4π6≤θ≤5π6,bbmax=0.076θ4-0.54θ3+1.1θ2-0.23θ-0.0056bbmax=0.076θ4-0.54θ3+1.1θ2-0.2π6≤θ≤5π6.The relative crack area (i.e., ratio of the partially open crack area to the fully open crack area) with respect to the angle of rotation is used as shown in Figure14. The results are similar when the crack closes. The relative crack area of simulation results is in good agreement with the crack closure straight line model [33] between rotation angles π/3 and 5π/6 and with the bilinear model in the range from 0 to π, but there are some relevant differences with respect to the crack closure perpendicular line used by Darpe et al. [12]. The area moment of inertia about the rotation axis can be approximated. Finally, the stiffness of the cracked shaft can be evaluated (stiffness is found to vary linearly with the area moment of inertia). The normalised stiffness (ratio between stiffness of the cracked shaft and stiffness of the uncracked shaft) is presented as shown in Figure 15. It is shown that the normalized stiffness between breathing mechanism results and the bilinear crack closure model is very close, which validates the accuracy of bilinear model.Figure 14
Comparison of relative crack area between breathing mechanism results, straight line, perpendicular, and bilinear crack closure model for a shallow cracka/d=0.1.Figure 15
Comparison of normalized stiffness between breathing mechanism results, straight line, perpendicular, and bilinear crack closure model for a shallow cracka/d=0.1.The fracture mechanics analysis presupposes the existence of an infinitely sharp crack leading to the singular crack tip fields. However, in real materials neither the sharpness of the crack nor the stress levels near the crack tip region can be infinite [27]. As an alternative approach to this singularity-driven fracture approach, the concept of CZM is proposed. CZM has evolved as a preferred method to analyze fracture problems not only because it avoids the singularity, but also because it can be easily implemented in a numerical method of analysis as in FE method. CZM presupposes the presence of a fracture process zone where the energy is transferred from external work in the crack regions. In this study, the external work flows as recoverable elastic strain energy and cohesive energy are examined, the latter encompassing the work and other energy consuming mechanisms within the fracture process zone. Process zone is defined as the region within the separating surfaces where the surface traction values are nonzero. This also implies that processes occurring within the process zone are accounted for only through the traction-displacement relations. The traction-separation relations for crack interfaces are demonstrated such that with increasing interfacial separation, the traction across the interface reaches a maximum; it means the elastic limit, where plastic dissipation in the surrounding (bounding) material is not accounted for in the process zone.


## 4.1. Breathing Mechanism Results

Figure9 displays the result of breathing mechanism during half-revolution of a rotating shaft (from closed to open crack) with relative crack depth  a/d=0.1. As can be seen the crack will start to appear and begin to open at π/6. The crack opens more slowly at the beginning, but increases its opening at π/3. At 5π/6 it is already completely open. Similarly, the crack closes again at 7π/6 and increases its closing at 3π/2. The crack is already completely closed at 11π/6. It is interesting to note that the angular positions where the crack opens again (from fully closed) are the same as where the crack begins to close.Figure 9
Breathing crack of the rotating shaft: relative crack deptha/d=0.1 at 600 rpm.Results of the relative crack depth versus the shaft rotation angle during breathing are shown in Figure10. Four functions are used to approximate the relative crack depth, namely, linear, quadratic, harmonic, and power functions. Table 1 gives the approximations for a relative crack depth a/d as a function of shaft rotation angle θ. It is shown that harmonic and quadratic functions have the minimal error and the nearest correlation with the simulation results. In this sense, the relative crack depth during crack opening should be understood not as linear increasing but as harmonic function or quadratic polynomial.Table 1
Curve fitting of the simulation results fora/d=0.1, θ in [rad].
Model
Mathematical model
SSE
R-square
Linear function
a
/
d = 0.039θ + 0.001
0.000 580
0.951 6
Quadratic polynomial
a
/
d = −0.012θ2 + 0.076θ − 0.020
0.000 011
0.999 0
Sinus function
a
/
d = 0.010sin (0.593θ − 0.145)
0.000 006
0.999 5
Power function
a
/
d = 0.123θ0.373 − 0.077
0.000 099
0.991 7Curve fitting of breathing mechanism results.
(a)
(b)
(c)
(d)


## 4.2. Discussion

Due to the presence of gravity, the upper portion of the cracked rotor at the beginning of a rotation is under compression and the crack is closed. As the rotor continues to rotate and the gravity direction is constant, the upper part now comes in the lower tensile region causing the crack to open. Moreover, the accuracy of the results is modelled by a bilinear model as shown in Figure11. The crack begins opening at shaft rotation angle θ0 and at crack front angle α that defines the angle between two crack front lines. The shaft rotates with shaft rotation angle θ and the crack front angle α increases until maximum crack width bmax is reached. The crack front angle α and the crack width b are nonlinearly increased during rotation from θ0 to θmax (Figures 12 and 13). One can obtain the crack front angle α and the crack width b as functions of the shaft rotation angle θ, respectively, by
(9)α={0,θ≤θ0crackstillclosed,-0.018θ3+0.016θ2+0.35θ+2.4,θ0≤θ≤θmaxcrackopens,0,θ≥θmaxcrack opens completely,bbmax={0,θ≤θ0crack still closed,0.076θ4-0.54θ3+1.1θ2-0.23θ-0.0056,θ0≤θ≤θmaxcrack opens,0,θ≥θmaxcrack opens completely.Figure 11
Bilinear crack closure model.Figure 12
Crack front angleα during rotation of shaft.Figure 13
Normalized crack widthb/bmax during rotation of shaft.For the relative crack deptha/d=0.1, the crack begins to open at  θ0=π/6  and  α=5π/6, while the crack opens completely at  θmax=5π/6  and  α=π. Therefore, the crack front angle and the crack width during crack opening are
(10)α=-0.018θ3+0.016θ2+0.35θ+2.4π6≤θ≤5π6,bbmax=0.076θ4-0.54θ3+1.1θ2-0.23θ-0.0056bbmax=0.076θ4-0.54θ3+1.1θ2-0.2π6≤θ≤5π6.The relative crack area (i.e., ratio of the partially open crack area to the fully open crack area) with respect to the angle of rotation is used as shown in Figure14. The results are similar when the crack closes. The relative crack area of simulation results is in good agreement with the crack closure straight line model [33] between rotation angles π/3 and 5π/6 and with the bilinear model in the range from 0 to π, but there are some relevant differences with respect to the crack closure perpendicular line used by Darpe et al. [12]. The area moment of inertia about the rotation axis can be approximated. Finally, the stiffness of the cracked shaft can be evaluated (stiffness is found to vary linearly with the area moment of inertia). The normalised stiffness (ratio between stiffness of the cracked shaft and stiffness of the uncracked shaft) is presented as shown in Figure 15. It is shown that the normalized stiffness between breathing mechanism results and the bilinear crack closure model is very close, which validates the accuracy of bilinear model.Figure 14
Comparison of relative crack area between breathing mechanism results, straight line, perpendicular, and bilinear crack closure model for a shallow cracka/d=0.1.Figure 15
Comparison of normalized stiffness between breathing mechanism results, straight line, perpendicular, and bilinear crack closure model for a shallow cracka/d=0.1.The fracture mechanics analysis presupposes the existence of an infinitely sharp crack leading to the singular crack tip fields. However, in real materials neither the sharpness of the crack nor the stress levels near the crack tip region can be infinite [27]. As an alternative approach to this singularity-driven fracture approach, the concept of CZM is proposed. CZM has evolved as a preferred method to analyze fracture problems not only because it avoids the singularity, but also because it can be easily implemented in a numerical method of analysis as in FE method. CZM presupposes the presence of a fracture process zone where the energy is transferred from external work in the crack regions. In this study, the external work flows as recoverable elastic strain energy and cohesive energy are examined, the latter encompassing the work and other energy consuming mechanisms within the fracture process zone. Process zone is defined as the region within the separating surfaces where the surface traction values are nonzero. This also implies that processes occurring within the process zone are accounted for only through the traction-displacement relations. The traction-separation relations for crack interfaces are demonstrated such that with increasing interfacial separation, the traction across the interface reaches a maximum; it means the elastic limit, where plastic dissipation in the surrounding (bounding) material is not accounted for in the process zone.


## 5. Conclusions

In this study, a simulation process of FE and MBS is employed. An elastic cracked shaft is modelled by FE which is transferred into MBS model. The vibration responses in the centroid of the shaft obtained from MBS have been transferred again into an FE model in order to observe the breathing mechanism and to study the stiffness losses due to breathing at different rotating speeds.The results have shown that the breathing mechanism is influenced by the vibration due to inertia forces and by relative crack depth. It is shown that the relative crack depth during crack opening should be understood not as linear increasing but as harmonic function or quadratic polynomial. It can be noted that as long as the relative crack depth is small, the bilinear crack closure model may be used. The main difference with respect to the straight line crack closure model is that the opening crack in the simulation is not constant at the beginning. The relative crack area of breathing mechanism results is in good agreement with the straight line crack closure model between rotation anglesπ/3 and π/6 and with the bilinear crack closure model in the range from 0 to π. Furthermore, the area moment of inertia about the rotation axis can be approximated. From this information, the stiffness losses of the shaft can be evaluated. It has been shown that the normalized stiffness between simulation and the bilinear crack closure model is very close, which validates the bilinear model.CZM are being increasingly used to simulate crack. Instead of an infinitely sharp crack envisaged in fracture mechanics, CZM presupposes the presence of a fracture process zone where the energy is transferred from external work in the crack regions. The traction-separation relations for crack interfaces are demonstrated such that with increasing interfacial separation, the traction across the interface reaches a maximum; it means the elastic limit, where plastic dissipation in the surrounding (bounding) material is not accounted for in the process zone.In light of the presented results and the conclusions, several subjects can be recommended for future research. An important point on which the knowledge could be improved is the prediction of crack propagation on cracked rotor and residual life estimation. The CZM can be easily implemented in FE model to analyze the dynamic behaviour of a cracked shaft. It is recommended to use CZM to study different types of cracks such as longitudinal and slant cracks and also multiple cracks. Some other parameters such as internal damping and thermal transients could be studied to obtain results for their effect on the breathing mechanism. Further analysis on crack morphology is extremely important to understand the dynamic behaviour of cracked shaft; this would include shallow and wide cracks.



---

*Source: 249035-2013-03-13.xml*