Türbin diskinin doğrusal olmayan gerilme gradyanı etkisi altındaolasılıksal çatlak ilerleme ömrü kestirimi
Probabilistic crack propagation of a turbine disc under the effectof nonlinear stress gradient
Yazar: ÖZGE TÜMERGİN
Danışman: PROF. DR. ATA MUĞAN
İstanbul Teknik Üniversitesi / Lisansüstü Eğitim Enstitüsü / Makine Mühendisliği Ana Bilim Dalı / Katı Cisimlerin Mekaniği Bilim Dalı
Derece: Yüksek Lisans, 2024
Konu: Makine Mühendisliği = Mechanical Engineering

Özet:
Gaz türbinli motorlarda katastrofik arızaları önlemek için tüm risk olasılıklarının
değerlendirilmesi gerekmektedir. Gaz türbinli motorlarda hasarın ana nedenlerinden
biri kırılma hasarıdır. Bu hasar ile ilgili tüm risklerin ortaya koyulabilmesi için
hesaplamalarda deterministik bir yaklaşım yerine stokastik bir yaklaşım kullanılması
önemlidir. Kırılma hasarı hesaplamalarında kullanılan ve parçanın ömrü üzerinde
doğrudan etkisi olan bazı parametreler, üretim veya test koşullarından kaynaklanan
belirsizlikler nedeniyle mutlak bir değerden ziyade istatistiksel bir dağılıma sahiptir.
Bu parametreler için mutlak sayısal değerler kullanıldığı durumlarda, tek bir çatlak
ilerleme ömrü sonucu elde edilecektir. Bu yaklaşım, aynı parametrenin farklı bir
olasılığa sahip diğer bir sayısal değerinden elde edilebilecek ömür sonucunu göz ardı
eder. Bu ömür sonucunun göz ardı edilmesi, sistemi katastrofik kırılma hasarı
riskiyle karşı karşıya bırakabilir.
Bu tez, olasılıksal çatlak ilerlemesi hesaplamalarını gerçekleştirmek için Matlab
kodu ile simüle edilen bir metodoloji sunmaktadır. Önerilen yaklaşımdan elde edilen
ömür sonuçları Nasgro çatlak ilerleme hesaplama programından elde edilen
sonuçlarla karşılaştırılmıştır. Daha sonra literatürden elde edilen deterministik ve
stokastik parametre değerleri yazılan olasılıksal çatlak ilerleme hesaplama koduna
dahil edilmiş ve elde edilen ömür dağılımı sonuçları, literatürde önerilen tasarım
ömrü sonuçları ile karşılaştırılmıştır. Bu bölümde önerilen yaklaşımda kullanılan
metodoloji özetlenecektir. İlk olarak, deterministik çatlak ilerlemesi hesaplama
metodolojisi anlatılacak, ardından olasılıksal yöntemin nasıl dahil edildiği
açıklanacaktır.
Bu tez çalışmasında Inconel 718 malzemeden imal edilen türbin motorunun şaft ve
disk parçaları ele alınmaktadır. Bu parçalar için olasılıksal çatlak ilerleme
hesaplamaları Matlab kodu kullanılarak gerçekleştirilmiştir. Çatlak ilerleme
hesaplama parametrelerinden biri olan gerilme değerinin elde edilebilmesi için türbin
modelinin sonlu eleman modeli Ansys programı kullanılarak oluşturulmuştur. Test
koşullarına uygun olacak şekilde, sıcaklık, basınç, tork ve açısal hız değerleri modele
sınır koşulları olarak tanımlanmıştır. Gerilme analizi, çatlak ilerlemesine neden
olabilecek en yüksek gerilme değerlerinin türbin diskinin civata deliği (bolthole)
kısmında gözlendiğini ortaya koymuştur. Çatlak ilerleme hesaplamaları bu kritik
bölge için gerçekleştirilmiştir.
Hasar kriterleri, kritik gerilme yoğunluğu faktörünün kırılma tokluğu değerini ve net
kesit gerilmesinin akma ve çekme dayanımı değerlerinin ortalamasını aşması olarak
belirlenmiştir. Bu kriterler gerçekleştiği andaki çevrim sayısı, parçanın çatlak
ilerleme ömrünü vermektedir. Çatlak ilerleme hesaplamalarında, hasar kriterlerini
xxiv
oluşturan gerilme yoğunluk faktörü ve net kesit gerilmesi değerlerinin hesaplanması
gerekmektedir.
Gerilme yoğunluğu faktörü, çatlağın tipine, geometrik oranlara ve parçanın gerilme
geçmişine bağlı bir değerdir. Gerilme yoğunluk faktörü hesaplamalarında bu üç etki
de dahil edilmiştir. Bu etkilerin ortaya koyulabilmesi için öncelikle kritik bölge
belirlenmelidir. Kritik bölge olarak, türbin diskinin civata deliği bölgesi seçilmiştir.
Gerilme yoğunluk faktörü hesaplamalarına dahil edilmesi gereken ilk etki olan çatlak
tipinin seçimi için öncelikle kritik bölge olarak kabul edilen türbin diskinin civata
bölgesi, çatlak ilerleme hesaplamalarını kolaylaştırmak için delikli bir plaka olarak
modellenmiştir. Çatlak tipi olarak plaka üzerindeki delik etrafında çeyrek eliptik köşe
çatlağı seçilmiştir.
Uygun çatlak tipi seçildikten sonra, gerilme yoğunluğu faktörünün hesaplanmasında
bir sonraki adım ikinci etki olan geometrik oranların belirlenmesidir. Plaka
üzerindeki delik etrafında çeyrek eliptik köşe çatlağı çatlak tipinin üç temel
geometrik oranı vardır: kalınlık yönündeki çatlak derinliğinin genişlik yönündeki
çatlak uzunluğuna oranı, kalınlık yönündeki çatlak derinliğinin plaka kalınlığına
oranı ve delik yarıçapının kalınlığa oranı. Yalnızca geometrik oranlara bağlı olan
referans gerilme yoğunluğu faktörleri, bu üç oranın belirli değerleri için
tablolaştırılmıştır. Geometrik oranlar, çatlağın ilerlemesine bağlı olarak her
döngünün sonunda değişmektedir. Dolayısıyla her döngünün sonunda, referans
gerilme yoğunluk faktörü değerleri de değişmelidir. Ana çatlak ilerleme koduna bir
alt kod eklenerek her döngü için yeni elde edilen geometrik oran için referans
gerilme yoğunluk faktörü değerleri hazırlanan tablodaki değerlerden interpole
edilerek hesaplanmıştır.
Gerilme yoğunluğu faktörünün hesaplanmasında ikinci adım tamamlandıktan sonra
sıra üçüncü etki olan gerilme geçmişine gelir. Çalışma alanı deliğin etrafında
olduğundan, gerilme konsantrasyonu etkisi de dikkate alınmalıdır. Gerilme
konsantrasyonu etkisi nedeniyle çatlağa uygulanan gerilme, deliğe en yakın konumda
en yüksek değerine ulaşır ve delikten uzaklaştıkça gerilme konsantrasyonu faktörüne
uygun olarak azalır. Bu durumu doğru bir şekilde modellemek için gerilme gradyanı
etkisinin hesaplamalara dahil edilmesi gerekir. Bu amaçla, gerilme yoğunluğu
faktörü değeri ağırlık fonksiyonu yöntemi kullanılarak hesaplanmıştır. Bu yöntem,
geometrik oranlara bağlı olarak değişen referans gerilme yoğunluğu değerlerini,
çatlağın plaka üzerindeki konumuna bağlı olarak değişen gerilme gradyanı ifadesi ile
integral kullanarak birleştirmektedir.
Çatlak ilerleme ömrünün belirlenmesinde kullanılan hasar kriterlerden biri olan
gerilme yoğunluk faktörünün hesaplanmasının ardından ikinci hasar kriteri olan net
kesit gerilmesinin hesaplanmasına geçilmiştir. Net kesit gerilmesi, kırılma hasarına
uğramış bir yapıda yerel akma ve ardından plastik çökme meydana gelip
gelmeyeceğini değerlendirmek için elastisite teorisi kullanılarak belirlenir.
Bu tez için yazılan Matlab kodu döngü bazında çalışmaktadır. Çatlak ilerlemesini ve
ilk çevrimi başlatmak için, hem genişlik hem de kalınlık yönlerindeki başlangıç
çatlak uzunluğu değerleri ve plakanın geometrik boyutları girilmelidir. Bu değerlere
dayanarak, plaka üzerindeki delik etrafındaki çeyrek eliptik köşe çatlağı için gerekli
geometrik oranlar hesaplanır. Gerilme değeri, çatlak ucunun mevcut konumuna bağlı
olan gerilme gradyanı ile belirlenir. Daha sonra, ağırlık fonksiyonu yöntemi
kullanılarak, gerilme yoğunluk faktörü değeri hem genişlik yönündeki çatlak hem de
derinlik yönündeki çatlak için ayrı ayrı hesaplanır. Bunun nedeni, çeyrek eliptik köşe
xxv
çatlağı durumunda, çatlağın her iki yönde de ilerlemesidir. Son olarak, net kesit
gerilme değeri de hesaplanır. Gerilme yoğunluğu faktörü ve net kesit gerilme değeri
hasar kriteri olan sınır değerler ile karşılaştırılır. Bu değerlerden herhangi biri kendi
sınır değerlerini aşarsa kod hesaplamaları durdurur. Eğer bu sınır aşılmazsa kodda
ikinci döngü başlar. Nasgro çatlak ilerleme denklemi, bir döngü için çatlak
uzunluğunun büyümesini hesaplar. Elde edilen uzama, başlangıçtaki çatlak uzunluğu
değerine eklenir. Bu yeni değer daha sonra bir sonraki döngü için yeni geometrik
oranları ve referans gerilme yoğunluğu faktörü değerlerini hesaplamak için kullanılır.
Çatlak büyüdükçe, plaka üzerindeki çatlak ucu konumu değişir ve kesit alanı azalır.
Bu da yeni çatlak uzunluğu için net kesit gerilimi değerinin yeniden hesaplanmasına
yol açar. Bu işlem, gerilme yoğunluğu faktörü değeri veya net kesit gerilme değeri
sınır değerlerini aşana kadar tekrarlanır. Bu değerlerden herhangi biri aşılırsa
hesaplamalar durdurulur ve mevcut çevrim sayısı çatlak ilerleme ömrü olarak kabul
edilir.
Bu noktaya kadar, deterministik çatlak ilerlemesi hesaplama metodolojisi
açıklanmıştır. Olasılıksal hesaplamalara geçilmeden önce sabit ve lineer olmayan
gerilme durumları için Nasgro programı ile yazılan kod sonuçlarının karşılaştırma
çalışması yapılmıştır.
Özetin önceki bölümünde, çatlak ilerleme hesaplamalarında kullanılan bazı
parametrelerin olasılıksal bir dağılıma sahip olması nedeniyle çatlak ilerleme
ömrünün de olasılıksal bir dağılıma sahip olması gerektiği belirtilmiştir.
Deterministik çatlak ilerleme algoritmasını olasılıksal çatlak ilerleme algoritmasına
çevirmek için Monte Carlo Simülasyonu yöntemi kullanılmıştır. Monte Carlo
simülasyonu, belirlenen iterasyon sayısında, rastgele değişkenler ile tekrarlı
hesaplamalar yapılması yoluyla, olasılıksal dağılıma sahip parametrelerin
deterministik hesaplamalarda kullanılmasını sağlayan bir yöntemdir. Literatürde,
10000 iterasyon yapılması önerilmektedir.
Bu tez, tahribatsız testler (NDT) ile belirlenen ilk çatlak boyu parametresini ve
malzeme testleri ile belirlenen kırılma tokluğu parametresini olasılıksal dağılıma
sahip parametreler olarak kabul etmektedir.
İlk çatlak boyu parametresinin literatürden kümülatif dağılım fonksiyonu grafiği elde
edilmiştir. Elde edilen grafikteki veriler kullanılarak ilk çatlak boyu dağılımına ait
ortalama ve standart sapma değerleri elde edilmiş daha sonra monte carlo
simülasyonlarında kullanılmak üzere 10000 adet rastgele değişken üretilmiştir.
Diğer olasılıksal dağılıma sahip olan parametre olan kırılma tokluğu için, test
sonuçları literatürden alınarak, ortalama ve standart sapma değerleri elde edilmiş ve
10000 rastgele değişken üretilmiştir.
Elde edilen rastgele değişkenlerin Monte Carlo simülasyonları kullanılarak
deterministik çatlak ilerleme hesaplamalarına dahil edilmesi sonucunda 10000 adet
çatlak ilerleme ömrü değeri elde edilmiştir. Ömür testinin sonuçları olasılıksal bir
dağılım grafiğine dönüştürülerek türbin diskinin çatlak ilerleme ömrünün risk
dağılımı ortaya çıkarılmıştır. Elde edilen risk dağılımı, literatürde bulunan güvenli
tasarım ömür hedefi ile karşılaştırılmıştır.

Abstract:
To prevent catastrophic failures in gas turbine engines, it is necessary to evaluate all
risk probabilities. The main cause of damage to gas turbine engines is fracture
failure, and it is important to calculate this using a stochastic approach rather than a
deterministic one to reveal all the risks. Calculations for fracture failure should be
carried out using the linear elastic fracture mechanics approach for cases that do not
include nonlinearity. By using this approach in conjunction with crack propagation
equations, we can obtain the crack propagation life for the part under examination.
Some of the parameters used in these calculations have a statistical distribution rather
than an absolute value due to uncertainties arising from production or test conditions.
These parameters have a direct impact on the life of the part. If absolute numerical
values are used to consider these parameters without taking into account their
statistical distributions, a single crack propagation life result will be obtained. This
approach also ignores the highly probable life result that would be obtained from
another numerical value of the same parameter with high probability. By ignoring
this life result, the outcome could put the system at risk of catastrophic fracture
damage.
This thesis presents a methodology for performing probabilistic crack propagation
calculations, simulated by a Matlab code. The life results obtained from the proposed
approach were compared with those from the Nasgro crack growth calculation
program. The methodology used in the proposed approach is summarised in this
section. First, the deterministic crack propagation calculation methodology is
explained, followed by an explanation of how the probabilistic method is
incorporated.
This thesis study discusses the shaft and disc parts of the turbine engine made of
Inconel 718 material. Probabilistic crack propagation calculations for these parts
were performed using Matlab code. To obtain the stress value, which is one of the
inputs, the finite element model of the components was analyzed using Ansys
software. Boundary conditions were applied to the model using temperature,
pressure, torque, and angular velocity values to match the test conditions. Linear
elastic analysis was performed since stress and temperature values were below the
limit values of the mechanisms (plasticity, creep) that could generate nonlinearity.
The stress analysis revealed that the highest stress values, which could cause crack
propagation, were observed in the bolthole part of the turbine disc. Crack
propagation calculations were carried out for the critical region due to its
susceptibility to fracture damage.
xxviii
The calculations for crack propagation utilized linear elastic fracture mechanics and
Nasgro crack growth equations. The damage criteria are based on exceeding the
critical stress intensity factor and flow stress.
The stress intensity factor measures the stress intensity at the crack tip. Fracture
damage occurs in the part when this value reaches the fracture toughness value,
which is a material property. The stress intensity factor is a value that depends on the
type of crack, geometric ratios, and stress history of the part. To obtain the stress
intensity factor value, it is necessary to include these three effects in the calculations.
The turbine disc, which is considered the critical region, was modelled as a plate
with a hole to facilitate crack propagation calculations. A quarter-elliptic corner
crack at the bolt hole was chosen as the crack type in accordance with the maximum
stress location.
Once the appropriate crack type has been selected, the next step in calculating the
stress intensity factor is to determine the geometric proportions. The bolthole crack
type has three fundamental geometric ratios: the ratio of crack length in the thickness
direction to the crack length in the width direction, the ratio of crack length in the
thickness direction to the thickness, and the ratio of bolthole radius to thickness.
Reference stress intensity factors, which depend solely on geometric proportions, are
tabulated for certain values of these three ratios. These factors can be used for all
cracks that have these geometric ratios. The geometric ratios change at the end of
each cycle due to the propagation of the crack. A sub-code has been added to the
main crack propagation code to calculate the reference stress intensity factor values
for the newly obtained geometric ratio by interpolating them from the values in the
prepared table. The stress history is then taken into account, considering the effect of
geometric proportions.
Once the second step in the calculation of the stress intensity factor has been
completed, it is time for the third step, the stress history. As the working area is
around the hole, the stress concentration effect must also be considered. Due to the
stress concentration effect, the stress applied to the crack reaches its highest value at
the position closest to the hole and decreases in accordance with the stress
concentration factor as it moves away from the hole. Therefore, the stress value that
causes crack propagation varies along the plate. To accurately model this situation,
the stress gradient effect must be included in the calculations. To this end, the stress
intensity factor value was calculated using the weight function method. This method
combines the reference stress intensity values, which vary with geometry, with the
stress gradient expression, which varies with the location of the crack on the plate,
using the integral.
After calculating the stress intensity factor, which is one of the main criteria for
determining the crack propagation life, the calculation of the net section stress, the
second criterion, was started. The net section stress is determined using strength of
materials theory to assess whether local yielding and subsequent plastic collapse will
occur in a fractured structure.
If the stress intensity factor exceeds the fracture toughness value or the net section
stress exceeds the flow stress, which is assumed to be the average value of the yield
and ultimate stress, the computation of the crack growth terminates.
The Matlab code written for this thesis works on a cycle-by-cycle basis. To initiate
crack growth and the first cycle, the initial crack length values in both the width and
xxix
thickness directions, as well as the geometric dimensions of the plate, must be
entered. Based on these values, the necessary geometric ratios for the quarter-elliptic
corner crack at the hole crack type are calculated. As previously stated, reference
stress intensity factor values are obtained through the interpolation method based on
these ratios. The stress value is determined by the stress gradient profile, which
depends on the current position of the crack tip. Then, using the weight function
method, the stress intensity factor value is calculated separately for both the crack in
the width direction and the crack in the depth direction. This is because in the case of
a quarter elliptical corner crack, the crack propagates in both directions. Therefore,
the crack propagation should be calculated in both directions. Finally, the net section
stress value is also calculated. The stress intensity factor is compared to the fracture
toughness value, and the net section stress value is compared to the flow stress value.
If either of these values exceed their respective limits, the code will stop. If the limits
are not exceeded, a final check remains before the crack propagation phase can
begin. The crack must exceed a certain stress intensity factor value to start
propagating. If it does not exceed this value, the crack remains as it is and does not
progress. This value is called the threshold stress intensity factor. If the stress
intensity factor exceeds this value, the crack will start to propagate.
The Nasgro crack growth equation is used for this stage, with fracture toughness and
stress intensity factor values as parameters. The Nasgro crack propagation equation
calculates the growth of a crack length for one cycle. The resulting elongation is
added to the initial crack length value. These calculations are made for both
directions of the crack tip. This new value is then used to calculate new geometric
ratios and reference stress intensity factor values for the next cycle. As the crack
grows, the crack tip position on the plate changes, and the cross-sectional area
decreases. This leads to a recalculation of the net section stress value for the new
crack length. This process is repeated until the stress intensity factor value exceeds
the fracture toughness or the net section stress value exceeds the flow stress value.
Calculations will be stopped, and the current number of cycles will be considered as
the crack growth life if any of these values are exceeded.
Up to this point, the methodology for calculating deterministic crack propagation has
been explained. The crack propagation results produced by the written deterministic
crack propagation code, both under constant stress and under a non-linner stress
gradient, were compared with the crack propagation results of the Nasgro program.
The previous section of the summary noted that crack propagation life should have a
probabilistic distribution due to certain parameters used in crack propagation
calculations having a probabilistic distribution. This thesis considers the initial crack
length parameter determined by non-destructive tests (NDT) and the fracture
toughness parameter determined by material tests as parameters with a probabilistic
distribution.
Initial crack length measurements are performed by NDT. The literature review has
indicated the NDT methods of initial crack depth testing for turbine disc life
estimation and found that the least risky method in terms of fracture risk is the Eddy
current method. Based on this information, high-resolution Eddy current test results
were obtained from the literature for the bolthole section of a turbine disc made of
nickel-based superalloy material with natural fatigue cracks. The source presents the
test results for the initial crack length as a cumulative distribution function (CDF)
graph. The probability that X takes a value less than or equal to x is known in
xxx
probability theory and statistics as the CDF. Probability Density Function (PDF) was
obtained by taking the derivative of the CDF. A subcode was developed to enable the
use of the CDF data in the crack propagation code. The code initially calculates PDF
values by numerically differentiating the CDF data. A probability density function
(PDF) is a function that provides the relative likelihood of a random variable's value
being equal to a given sample in the sample space. It describes the set of possible
values that the variable could take. Thus, it was determined which values the initial
crack length parameter can take. The code generates a set number of initial crack
length values based on the given probabilities. It then calculates the mean and
standard deviation of these values. Finally, it generates 10.000 random initial crack
length values using the calculated values.
After obtaining the initial crack length values, it is time to obtain the fracture
toughness values, which are considered as a probabilistic parameter within the scope
of this thesis. Measurement results regarding the fracture toughness parameter of
Inconel 718 material were obtained from the literature. The mean and standard
deviation of these results were calculated. Subsequently, 10.000 random variables
were generated based on the calculated mean and standard deviation values.
Monte Carlo simulations were then used to model the probabilistic crack propagation
algorithm. Monte Carlo simulation is a numerical algorithm used to obtain results
through repeated random samples. The deterministic section mentions a cyclic
algorithm that was run 10.000 times using different previously produced random
values of initial crack length and fracture toughness parameters. This study resulted
in 10.000 crack propagation life values. The results of the life test were transformed
into a probabilistic distribution chart, revealing the risk distribution of the turbine
disc's crack propagation life.
The obtained life distribution results were compared with the design life results
suggested in the literature.