Anabilim Dalı: Makine Mühendisliği Programı: Katı Cisimlerin Mekaniği
ĠSTANBUL TEKNĠK ÜNĠVERSĠTESĠ FEN BĠLĠMLERĠ ENSTĠTÜSÜ
ÖZELLĠKLERĠ FONKSĠYONEL OLARAK DEĞĠġEN MALZEMELERDE GERĠLME DALGALARININ ĠNCELENMESĠ
YÜKSEK LĠSANS TEZĠ
Müh. Ġlyas TOPRAK
Tez DanıĢmanı: Doç.Dr. Erol ġENOCAK
ÖNSÖZ
Bu tez çalışmasındaki yardımlarından dolayı Doç. Dr. Erol Şenocak ve Hakan Tanrıöver’e, tez yazım süresince bilgisayarını kullanmama izin veren Sami Çıray, Ozan Durmazoğlu ve Doğan Akcan’a, şekiller konusunda yardım eden Emre Uzunoğulları ve İdris Tasak’a teşekkür ederim.
İÇİNDEKİLER ÖNSÖZ KISALTMALAR TABLO LİSTESİ ŞEKİL LİSTESİ ÖZET SUMMARY 1. GİRİŞ ii v vı vıı x xı 1 2. TEMEL DENKLEMLER
2.1 Sonlu Elemanlar Denklemlerinin Çıkarılışı 2.2 Eksenel Simetri Denklemleri
2.3 Fonksiyonel Olarak Değişen Malzemeler
3 3 6 7 3. SONLU ELEMANLARIN PROGRAMLAMA YÖNTEMİ
3.1 Koordinat Dönüşümleri
3.2 Entegrallerin Sayısal Olarak Hesaplanması
9 9 9 4.ZAMAN ENTEGRASYONLARI 4.1 Açık Yöntem 4.2 Kapalı Yöntem
4.3 Açık Zaman-Entegrasyon prosedürü 4.3.1 Merkezi Farklar Yöntemi
4.3.2 Merkezi Farklar Yöntemi Sonucunda Ortaya Çıkan Denklemin Çözümü
4.4 Kapalı Zaman-Entegrasyon Prosedürü 4.4.1 Newmark Yöntemi
4.4.2 Newmark Yötemi Sonucunda Ortaya Çıkan Denklemin Çözümü 11 11 12 12 12 14 15 15 16 5.SAYISAL ÖRNEKLER
5.1 Çubuğun Rijid Duvara Çarpması
5.2 Silindirik Çubuk Üzerine Gelen Anlık Bir Kuvvet Sonucu Oluşan Yer Değiştirmeler
5.3 Elastisesi Değişen Silindirik Plağın Statik Durumdayken Analitik Çözümleriyle Yazılan Programın Sonuçlarının Karşılaştırılması
17 17 22 25
5.4 Elastisesi Değişen Dairesel Bir Plağın Üzerine Gelen Basınç Etkisiyle Plak Üzerindeki Bir Noktanın Yer değiştirmesinin Zamana Bağlı Bulunması
5.5 Elastisitesi Değişen Silindirik Plağın Üzerine Gelen Anlık Bir Basınç Durumundaki Halinin İncelenmesi
28 33 6.SONUÇLAR VE TARTIŞMA 35 KAYNAKLAR EKLER 24 ÖZGEÇMİŞ 36 39 47
ISALTMALAR K :Katılık Matrisi M :Kütle Matrisi C :Sönüm Matrisi :Poisson oranı :Yoğunluk W :Ağırlık katsayısı w :Çökme Miktarı :Gerilme :Gerinme
TABLO LİSTESİ
Sayfa No
Tablo 3.1 n değerleri için nip, i,j wiwj değerleri... 10
Tablo 5.1 Modellemede kullanılan dikdörtgen çubuğa ait olan özellikler... 18
Tablo 5.2 Bulunan yer değiştirme değerleri ve referans [4] den alınan değerlerin farkları... 19
Tablo 5.3 Modellemede ve analizde kullanılan silindirik çubuğa ait olan özellikler... 24
Tablo 5.4 Analiz zamanları(Cpu zamanı)... 24
Tablo 5.5 h/a=0.05 için orta noktanın yer değiştirmesi... 26
Tablo 5.6 h/a=0.15 için orta noktanın yer değiştirmesi... 27
Tablo 5.7 h/a=0.2 için orta noktanın yer değiştirmesi... 27
ŞEKİL LİSTESİ
Sayfa No
Şekil 2.1 Genel üç boyutlu kütle... 3
Şekil 5.1 Düzlem gerinme halindeki çubuğun duvara çarpması... 17
Şekil 5.2 Çubuğun çözüm ait ve düğüm noktaları numaraları... 18
Şekil 5.3 42. dügüm noktasının yer değiştirmedi [4]... 19
Şekil 5.4 42. dügüm noktasının yer değiştirmedi ... 20
Şekil 5.5 Merkezi faklar yöntemi kullanılarak çubuk üzerinde oluşan 1sn sonraki yer değiştirmeler... 21
Şekil 5.6 Newmark yöntemi kullanılarak çubuk üzerinde oluşan 1sn sonraki yer değiştirmeler... 21
Şekil 5.7 Silindirik bir çubuğun üzerine gelen anlık yükler... 22
Şekil 5.8 Çubuğun yarısının çözüm ağı ve incelenen düğüm noktasınınn yeri... 23
Şekil 5.9 Uygulanan kuvvet ve geometrinin şekli... 25
Şekil 5.10 Malzemenin kesit boyunca, elastisitesinin n ile değişimi... 26
Şekil 5.11 Denklem (5.1)’in değişik “n” lerde, eleman sayısına göre yakınsaması... 28
Şekil 5.12 Analitik olarak çözülen geometri modeli... 31
Şekil 5.13 Yapının çözüm ağı ve incelenen 3. düğüm noktası ... 31
Şekil 5.14 Yapının merkezindeki yer değiştirmeleri... 32
Şekil 5.15 Newmark ve Merkezi Fark yöntemi Newmark ve Merkezi Fark yöntemi kullanılması sonucu yapının merkezinde oluşan yer değiştirmeler... 32
Şekil 5.16 Plağın eksen boyunca malzeme değişimi... 34
Şekil 5.17 Plak geometrisi... 34
Şekil A.1 Newmark yöntemi ile t=0.008 milisaniyedeki çubuğun üzerindeki eksenel yer değiştirmeler... 38
Şekil A.2 Merkezi farklar yöntemi ile t=0.008 milisaniyedeki çubuğun üzerindeki eksenel yer değiştirmeler... 38
Şekil A.3 Newmark metodu ile t=0.008 milisaniyedeki çubuğun üzerindeki radyal yer değiştirmeler... 39 Şekil A.4 Merkezi farklar metodu ile t=0.008 milisaniyedeki çubuğun
üzerindeki radyal yer değiştirmeler... 39 Şekil A.5 z=25.4 mm ve z=50.8 mm deki yerdeğiştirmeleri göstermektedir
[6]... 40 Şekil A.6 6 Merkezi Farklar Yöntemi kullanılarak değişik zaman
adımlarında ve değişik noktalardaki radyal yerdeğiştirmeler. ... 40 Şekil A.7 Newmark kullanılarak değişik zaman adımlarında ve değişik
noktalardaki radyal yerdeğiştirmeler... 41 ŞekilA.8 Newmark Yöntemi kullanılarak, 1e-4 zaman aralığında, 50x10
eleman sayısında, r=10 mm ve z=254 mm noktasındaki düğümün eksenel yer değiştirmesi... 41 Şekil A.9 Merkezi Farklar Yöntemi kullanılarak, 1e-4 zaman aralığında,
50x10 eleman sayısında, r=10 mm ve z=254 mm noktasındaki düğümün radyal yer değiştirmesi... 42 Şekil A.10 Merkezi Farklar Yöntemi kullanılarak, 1e-4 zamanında, 50x10
eleman sayısında, r=10 mm ve z=254 mm noktasındaki düğümün eksenel yer değiştirmesi... 42 Şekil A.11 Merkezi Farklar Yöntemi kullanılarak, 1e-4 zamanında, 50x10
eleman sayısında, r=10 mm ve z=254 mm noktasındaki düğümün eksenel yer değiştirmesi... 43 Şekil A.12 r=10 mm ve z=24.5 mm deki düğüm noktasının, Newmark
Farklar Yöntemi kullanılarak 0.1 milisaniye zaman aralığındaki radyal yer değiştirmesi... 43 Şekil A.13 r=10 mm ve z=24.5 mm deki düğüm noktasının, Newmark
Farklar Yöntemi kullanılarak 0.1 milisaniye zaman aralığındaki eksenel yer değiştirmesi... 44 Şekil A.14 r=10 mm ve z=24.5 mm deki düğüm noktasının, Merkezi Farklar
Yöntemi kullanılarak 0.1 milisaniye zaman aralığındaki radyal yer değiştirmesi... 44
Şekil A.15 r=10 mm ve z=24.5 mm deki düğüm noktasının, Merkezi Farklar Yöntemi kullanılarak 0.1 milisaniye zaman aralığındaki eksenel yer değiştirmesi... 45
ÖZET
Bu çalışmada mekanik özellikleri fonksiyonel olarak değişen malzemelerin özellikleri ve mekanik davranışları incelendi. Bu malzemelerin dinamik yükler altındaki deformasyon alanlarının belirlenmesi için gerekli denklemler ve sonlu elemanlar formülasyonları elde edildi. Açık (Merkezi Farklar) ve kapalı (Newmark) yöntemler kullanılarak elde edilmiş formülasyonlar çeşitli problemlere uyarlanarak literatürde bulunan sonuçlarla karşılaştırıldı. Bu çalışma boyunca nispeten basit sınır koşullarına sahip ve özellikleri fonksiyonel olarak değişen yapıların dinamik davranışları incelenmiştir. Yapının mukavemet özelliklerinin zamanla değişmediği ve sistemin tamamen doğrusal olduğu farz edilmiştir.
Bu çalışma içerisinde, silindirik bir çubuk ve plak için yapılan sonlu elemanlar analizlerinde Newmark Yöntemi, Merkezi Farklar Yöntemine göre daha iyi sonuçlar verdiği görülmüştür. Belirli bir hata sınırı için Newmark Yöntemi daha geniş zaman adımlarının kullanılmasını ve böylece sayısal çözümün çok daha kısa sürede sonuçlanmasını sağlamaktadır. Zaman adımının çok küçük seçilmek zorunda kalındığı koşullarda ise Merkezi Farklar Yönteminin çözüm süresi Newmark Yöntemine göre daha düşük çıkmıştır.
SUMMARY
In this thesis, properties and mechanical behaviors of functional graded materials were studied. Necessary equations and finite element formulations for the strain fields of these materials under dynamic loads were derived. Implicit (central difference) and explicit (Newmark) methods were utilized in the derived formulations and results were compared with literature. During the study, dynamic behaviors of functional graded materials with simple boundary conditions were analysed. Mechanical properties are assumed to be time independent and the system is assumed to be linear.
The study showed that the Newmark method used in the finite element analysis of cylindirical rods and plates gives better results compared to the central differences method. Under a certain degree of error, Newmark method allows larger time steps and hence decreases the amount of computation time. In cases where time step has to be much smaller, the central differences method gave faster results compared to the Newmark method.
1.GİRİŞ
Özellikleri yapı içerisinde noktadan noktaya değişen yapılara özellikleri fonksiyonel olarak değişen yapılar denir. Bu çalışmanın amacı bu tip yapıları incelemek ve dinamik yükler altındaki davranışlarını bulmaktır. Çözücü olarak açık ve kapalı yöntemler kullanılmıştır.
Son yüzyılda, özellikleri fonksiyonel olarak değişen malzemeler üzerine incelemelerin sayısı gözle görülür biçimde artmıştır. Bu tip malzemelerin incelenmelerinin başlıca amacı daha dayanıklı veya daha hafif yapılar elde etmektir. Bu tip yapılar askeri, havacılık ve otomotiv sanayilerinde gün geçtikçe daha geniş kullanım alanı bulacaktır. Örneğin, geleceğin tank zırhları bu tip malzemelerden yapılarak daha hafif bir yapıya sahip olacaklardır. Bu zırhları kullanan tanklar daha süratli ve saldırılara karşı daha dayanıklı hale geleceklerdir. İleri teknolojiye dayalı kompozit malzemeler kullanılarak aynı mukavemete sahip daha hafif uçak komponentleri imal edilecektir.
Kompozit malzemelerin üretiminden önce bilgisayar ortamında bu tip yapıların değişik yükler altında nasıl şekil değiştireceği modellenmelidir. Bu tip simülasyonların sonuçları kullanılarak belli bir yük dağılımına dayanabilecek ve mekanik özellikleri malzeme içerisinde değişen (fonksiyonel değişen) komponentler tasarlanabilir. Bu konudaki çalışmalar Berezovski, Engelbrecht ve Maugin [7], Erdoğan ve Chiu [1], Gazonas, Thamburaj ve Santare [8], Liu ve Han [10], Velasco ve Zarate [11] tarafından yayınlanan makaleler olarak gösterilebilir. Fonksiyonel malzemelerle ilgili yapılan yayın çalışmalarının son yıllarda gösterdiği artış ileride özellikle askeri amaçlı sanayilerde bu malzemenin geniş bir kullanım olanağı bulacağını göstermiştir.
Bu çalışmanın amacı, özellikleri fonksiyonel olarak değişen malzemeleri incelemek ve bu malzemelerin mekanik davranışları için analitik ve sayısal modeller geliştirmektir. Bu yapıların analitik yolla çözülmesi ancak basit sınır koşulları için yapılabilmektedir. Bu nedenle fonksiyonel değişen yapıların mekanik davranışları
ancak sonlu elemanlar simülasyonu ve benzeri sayısal yöntemlerle incelenebilmektedir. Sayısal yöntemlerin kullanılabilmesi için, fiziksel sistemin matematiksel modelinin zaman ve koordinat uzaylarında ayrıklaştırılması gerekmektedir. Zaman içerisinde ayrıklaştırma malzemenin dinamik davranışının modellenmesi için kritik bir aşamadır. Bu ayrıklaştırma açık (explicit) ve kapalı (implicit) olmak üzere iki yöntemle yapılmaktadır. Bu çalışmada açık yönteme ait Merkezi Farklar Yöntemi ve kapalı yönteme ait Newmark Yöntemi kullanılarak fonksiyonel değişen malzemelerin dinamik yükler altındaki davranışları incelenip sonuçlar karşılaştırılmıştır.
2. TEMEL DENKLEMLER
Bu bölümde genel denklemler çıkarılacak, daha sonra bu denklemler sonlu elemanlar yöntemine göre tekrar elde edileceklerdir. Sonlu elemanlar için çıkarılan denklemler, özel bir koşul olan eksenel simetri durumu için yeniden çıkarılacaktır.
2.1 Sonlu Elemanlar Denklemlerinin Çıkarılışı
Şekil 2.1‟de, denge halindeki genel üç boyutlu kütle görülmektedir. Kütle üzerine gelen kuvvetler yüzey kuvvetleri fS, kütle kuvvetleri fB, ve tekil kuvvetler Fi dir. Bu kuvvetler, kütle üzerine gelen tüm yükleri temsil etmektedirler.
Şekil 2.1 Genel üç boyutlu kütle
B Z B Y B X f f f B f , S Z S , Y S X f f f S f , i Z i Y i X F F F i F (2.1)
Kütle üzerine gelen kuvvetler (2.1)‟de verilmiştir.
U V W
T
U (2.2)
Yer değiştirmeler ise U vektörü olarak (2.2)‟de verilmiştir. Z,W X,U Y,V x, u z,w y,v Fzi Fxi Fyi i fzs , fzB fxs , fxB fys , fyB
zx yz xy zz yy xx
Tε
(2.3)
zx yz xy zz yy xx
Tτ
(2.4)Denklem (2.3) ve (2.4) ise sırasıyla gerinme ve gerilme vektörleridir.
Şekil 2.1 deki kütlenin denge durumunu yazmak için sistemin toplam potansiyel enerjisine bakacağız. Sistemde toplam iç ve dış enerji farkları sabit olmak zorundadır.
S i V V i T i S T S B T T 2 1 F U dS f U dV f U dV Πε
τ
(2.5)Bu ilkeyi kullanarak, denklem (2.5) kolaylıkla yazılabilir. Virtuel iş ilkesi doğrudan toplam potansiyel ilkesine uygulanabilir. Virtuel iş ilkesine göre, sistemde yapılacak küçük bir yer değiştirme sonucu, sistemde oluşacak toplam iç enerji, sistemde oluşan toplam dış enerjiye eşit olmalıdır.
ε
τ
D (2.6)
S i i T i S T S V B T V T δ δ δ δε
Dε
dV U f dV U f dS U F (2.7)Denklem (2.6) gerinme arasındaki ilişkiyi göstermektedir. D gerilme-gerinme matrisidir. Denklem (2.6), denklem (2.5) de yerine konur ve virtuel iş ilkesi uygulanırsa denklem (2.7) elde edilir.
ε
ε
; δUU; US US; δUi Ui (2.8) Denklem (2.7) deki ifadeler daha basit şekilde (2.8) deki gibi ifade edilebilir.
Şekil 2.1‟deki yapıyı, sonlu elemanlar yöntemi için gerçek yapıya yaklaşık, sonlu elemanlara ayrılmış ve düğüm noktalarından birleştirilmiş bir yapı olarak düşünebiliriz. Her eleman için yer değiştirmeler düğüm noktalarında hesaplanır ve her elemanın yer değiştirmesi, yer değiştirmenin bir fonksiyonu olarak yazılabilir. Bir “m“ elemanı için hesap yapacak olursak. Yer değiştirme
U N
u(m)(x,y,z) (m)(x,y,z)ˆ (2.9) denklem (2.9)‟daki gibi olur. N(m) burada yer değiştirme interpolasyon matrisidir. Uˆ
ise 3 boyutlu yer değiştirme vektörüdür. N sonlu eleman için Uˆ
1 1 1 2 2 2 N N N
T W V U .... W V U W V U ˆ U (2.10)
1 2 N
T U .... U U ˆ U (2.11)Her düğüm noktası için birim uzama yer değiştirme bağıntısı ise şu şekilde yazılabilir. U B
ε
(x,y,z) (x,y,z)ˆ (m) (m) (2.12)Denklemde birim uzama yer değiştirme matrisi, B(m), N(m) matrisinin türetilmiş halidir.
En genel haldeki hareket denklemi (2.1) deki gibi yazılabilir. Şimdi ise elemanların gerilme-gerinme ilişkisi gösterelim.
I m m m m) ( ) ( ) ( ) (
τ
τ
D ε (2.13)D(m) katılık matrisi,
τ
I(m) ise öngerilmedir. Her eleman için denklemler tekrarçıkarıldıktan sonra, denklemler (2.13), (2.12), (2.9) ve (2.8), denklem (2.7)‟de yerine yazılacak olursa denklem (2.14) elde edilir.
F dV B dS f N dV f N U U dV B D B Uτ
m V (m) I (m) T (m) m S (m) S (m) S (m) m V (m) B (m) T (m) T m V (m) (m) (m) T (m) T (m) (m) (m) (m) ˆ ˆ ˆ (2.14) KU=R R=RB+RS-RI+RC (2.15) (2.16) Denklem (2.14) denklem (2.15)‟deki gibi yazılabilir. Burada R kuvvet vektörü, K katılık matrisi, U ise yer değiştirme vektörüdür.
m V (m) (m) (m) T (m) (m) dV B D B K
m V (m) B (m) T (m) B (m) dV f N R
m S (m) S (m) S (m) S (m) dS f N R Rc=F (2.17) (2.18) (2.19) (2.20)Denklem (2.15) ve denklem (2.16)‟nın açık ifadeleri denklem (2.17), (2.18), (2.19) ve (2.20)‟de verilmiştir.
En genel haldeki hareket denklemi (2.21) deki gibi yazılabilir.
t (t) (t) (t) KU R U C U M (2.21)Denklemdeki K ve R matrislerinin açık hali daha önceden gösterilmiştir. C sönüm ve M kütle matrisleri ise D‟Albembert prensibi kullanılarak bulunacaktır. Denklem (2.20) ye D‟Albembert prensibini uygulayacak olursak
m V (m) (m) (m) (m) (m) B (m) T (m) B (m) dV U N κ U N ρ f N R (2.22)RB denklemi (2.22) halini alır.
m V (m) (m) T (m) (m) (M) dV N N ρ M
m V (m) (m) T (m) (m) (m) dV N N κ C (2.23) (2.24)Kütle ve sönüm matrisi ise denklem (2.23) ve (2.24) deki gibi gösterilebilir.
2.2 Eksenel Simetri Denklemleri
Bölüm 2.1‟de çıkarılmış olan denklemler bu bölümde eksenel simetri özelliği için tekrar gösterilecektir. Denklemler tekrar gösterilirken e‟inci eleman için gösterilecek ve sönümün olmadığı farz edilecektir.
dz dr r
dV (2.25)
Eksenel simetrik bir eleman için birim hacim, denklem (2.25) deki gibi ifade edilebilir. (2.25) denklemindeki r ve z silindirik koordinatları göstermektedir. (2.25) deki ifade denklem (2.17) ve denklem (2.23)‟de yerine yazılırsa
B DBrdrdz Ke Ωe T dz Nrdr ρN Me
Ωe T (2.26) (2.27) denklem (2.26) ve (2.27) elde edilir.
ε
Au w u r r z z r rz z r 0 1 0 0 (2.28)Eksenel simetrik bir durum için mühendislik gerinmesi (2.28) denkleminde gösterilmiştir. Burada u radyal ve w eksenel yöndeki yer değiştirmeyi göstermektedir. N şekil fonksiyonu matrisi olmak üzere B=AN şeklinde ifade edilebilir. 1 0 1 1 0 ) 1 ( 2 2 1 0 0 1 0 1 1 1 0 1 1 ) 2 1 )( 1 ( ) 1 ( E D (2.29)
Gerilme-gerinme matrisi ise (2.29) şeklindedir. Burada E elastisite modülünü, ve Poisson oranını gösterir.
2.3 Fonksiyonel Olarak Değişen Malzemeler
Özellikleri fonksiyonel olarak değişen malzemelerde elastisite modülü ve yoğunluk, r ve z nin birer fonksiyonudur. Bundan dolayı katılık ve kütle matrislerinde bulunan ve E entegral dışına çıkamazlar. Denklemler (2.26) ve (2.27) denklemler (2.30) ve (2.31) haline gelirler. Poisson oranındaki değişme ise ihmal edilmiştir. Gerilme-gerinme matrisinin yeni hali ise (2.32)‟de verilmiştir.
B D Brdrdz K T ( r ,z) e e dz Nrdr N ρ M T
e (r,z) e (2.30) (2.31) 1 0 ν 1 ν ν 1 ν 0 ν) 2(1 2ν 1 0 0 ν 1 ν 0 1 ν 1 ν 1 ν ν 0 ν 1 ν 1 2νν ν)(1 (1 ν) z)(1 (r, E D (2.32)
3. SONLU ELEMANLARIN PROGRAMLAMA YÖNTEMİ
Bu bölümde, ikinci bölümde oluşan denklemlerin sayısal olarak çözülebilmesi için gerekli olan dönüşümler gösterilmiştir.
3.1 Koordinat Dönüşümleri
Yazılacak olan programda kullanılan eleman tipi bilineer dikdörtgen elemandır. Elemanlar yerel koordinatlarda kullanıldıkları için bazı dönüşümlere ihtiyaç vardır. Bu dönüşümler zincir kuralı uygulanarak yapılır.
y x J y x y x (3.1)
Denklem (3.1) deki J, jakobiyan matrisini göstermektedir. Denklem (3.2), entegralin yerel koordinatlara dönüştürülmesini göstermektedir.
1 1 1 1 detJdd dxdy (3.2)3.2 Entegrallerin Sayısal Olarak Hesaplanması
Bir çok entegralin analitik olarak hesaplanması mümkün olmasına karşın, bilgisayar programları için analitik hesap uygun bir yöntem değildir. Bunun yerine bir çok sonlu elemanlar programında entegraller sayısal olarak hesaplanır.
(3.3)
nip i j i i n i n j j i j iw f W f w d d f 1 1 1 1 1 1 1 ) , ( ) , ( ) , ( entegralin sayısal yöntemlerle yaklaşık olarak nasıl hesaplandığı denklem (3.3)‟de verilmiştir. nip=n2
entegrasyon noktaların toplamı wi ve wj veya Wi ağırlık katsayısı j
i
, ise Gauss entegrasyon noktalarının koordinatlarıdır.
Tablo 3.1 n değerleri için nip, i,j wiwj değerleri.
n nip i,j wiwj 1 1 0 2 2 4 3 1 4 3 9 5 3 5/9
Tablo 3.1‟de n‟nin 1, 2 ve 3 değerleri için wiwj, i,j ve nip değerleri
gösterilmektedir [4]. Katılık matrisinin entegral şekli aşağıda verilmiştir.
Ke
eBTDBrdrdz
2 1 2 1 ) (det i j ij j iw J w BTDijBij (3.4)aynı şekilde kütle matrisi de yazılırsa
Me
eρNTNrdr
2 1 2 1 ) (det i j ij ij ij j iw J w NTDijN (3.5)4. ZAMAN ENTEGRASYONLARI
Bu bölümde açık ve kapalı yöntemlerin neler olduğu, hangi tip analizlerde kullanıldıkları ve formülasyonları verilecektir. Daha sonra kapalı ve açık yöntemler için türetilen denklemler hareket denkleminin içerisine yerleştirilecek ve oluşan denklemlerin nasıl çözüldüğü anlatılacaktır.
4.1 Açık Yöntem
Açık yöntemler, kısa süreli dinamik olaylar veya dalga yayılımı problemleri için çok uygundur. Buna örnek olarak patlamalar veya bir yere çok hızlı çarpma eylemleri gösterilebilir. Açık yöntemlerde tt deki sonuçları bulabilmek için t anındaki denge denklemlerinin kullanılması gerekmektedir. Her adım adım çözümde katılık matrisinin tekrar çözüme girmesine gerek yoktur ve diyagonal kütle matrisi kullanılırsa, ek bir depolamaya da gerek kalmaz. Adım başına hesaplama süresi ve depolama ihtiyacı kapalı yöntemlere göre çok daha azdır. Bunların yanında açık yöntemler kullanırken, iterasyonu kararlı halde tutmak için zaman aralıklarının çok küçük seçilmesi gerekmektedir. Açık yöntemlerde, zaman aralıkları genellikle çözüm hassasiyetinden çok daha önemlidir. Orta süreli dinamik problemler ele alındığında ise kullanılan bu kısıtlamalar, yaklaşımın etkinliğini sınırlandırmaktadır. Açık yöntemler tt zaman adımındaki çözümde bize t zamanı için yer değiştirme, hız ve ivmeyi hesaplamamıza izin verir. Bu sayede malzeme yapısı tam olarak değerlendirilerek, bir sonraki zaman adımı için yapılacak çözümdeki iç gerilmelerin hesaplanmasını sağlar. Bu da çok karışık malzeme modellerinin bile kolaylıkla çözüme dahil edilmesini sağlar Dokainish ve Subbaraj [4].
4.2 Kapalı Yöntem
İkinci yöntem ise kapalı zaman entegrasyon yöntemidir. Güçlü olduğu bölgeler kapalı yöntemlerin zayıf olduğu eylemsizlik, sismik problemleri ve karmaşık yapı geometrileri bulunan alanlardır. Kapalı yöntemde o anda hesaplanan yer değiştirme denklemi o andaki hız ve ivmeyi içerir. Bu yüzden tt adımındaki yer değiştirmenin bulunması, katılık matrisinin her zaman adımı için çözülmesini gerektirir. Kapalı yöntemlerin çoğu lineer çözümler için koşulsuz olarak kararlıdır ve kullanılacak en büyük zaman uzunluğu entegrasyon sürecinin kararlılığıyla değil çözümün kesinliğiyle ilgilidir. Kapalı yöntemler açık yöntemlere kıyasla her adımda daha fazla hesaplama gerektirse de, zaman aralığı, sadece kesinlik ölçütleri ile sınırlandırıldığı için çok daha geniş tutulabilir. Diğer yandan çoğu açık yöntemde zaman aralığı nümerik kararlılıkla sınırlı olduğundan, elde edilecek çözümler gerekenden daha kısa zaman dilimleri için hesaplandığından hesaplama maliyeti daha fazla olmaktadır.
Ne yazık ki açık ve kapalı yöntemlerin bir dezavantajı vardır. tt zamanındaki malzemedeki durum bir önceki adımdan tahmin edilmek zorundadır. Bu bazen çok güç veya imkansız olabilmektedir. Bu olay malzeme modelinin karmaşıklığı veya doğrusal olmaması nedeniyle meydana gelmektedir. Ayrıca her zaman problemden bağımsız olarak sayısal hatalar meydana gelmektedir
4.3 Açık Zaman-Entegrasyon Prosedürü
Yapısal analizlerde, açık çözüm tekniği t zaman adımı sonundaki yer değiştirme vektörünü verir. Yani t zaman adımı sonundaki yapısal yükler veya katılık matrisi önemli değildir. Bu iteratif hesaplama adım adım hesaplamalarda bize kolaylık sağlar. Açık zaman-entegrasyon yöntemlerinde hız ve ivmeyi bulabilmek için, bulunan yer değiştirmelerin farklı bir formülasyon içine konulması gerekmektedir.
4.3.1 Merkezi Farklar Yöntemi
Merkezi farklar yöntemi, yapısal dinamik problemlerinde çok yaygın olarak kullanılan bir açık yöntemdir. Merkezi farklar yöntemi diğer açık yöntemlere göre daha kararlı ve yüksek doğruluğa sahiptir. Ancak kötü bir yanı vardır. Çok küçük
zaman adımlarına ihtiyaç duyar. Merkezi farklar yöntemi aşağıdaki merkezi farklar formülasyonunu temel alır.
t
/
)
(
t Δt t 2 / Δt t
X
X
X
Δt / ) ( t t Δt 2 / Δt X X Xt t / ) ( t Δt/2 t Δt/2 t X X X Δt / ) ( t Δt/2 t Δt/2 t X X X 2 / ) ( t Δt t 2 / Δt t X X X ve 2 / ) ( t Δt 2 / Δt t X Xt X (4.1) (4.2) (4.3) (4.4) (4.5) (4.6) denklemlerde bulunan tt, t ve tt zamanları birbirlerini izleyen zaman seviyeleridir. Merkezi farklar yönteminin yapısını oluştururken, yukarıdaki denklemleri başlangıç olarak kullanacağız. (4.5) ve (4.6) numaralı denklemleri hızı bulma için (4.3) nolu denklemin içinde, ivmeyi bulmak için de (4.1) ve (4.2) nolu denklemleri (4.4) nolu denklemde t zaman adımında kullanacağız. XtΔt, X ve tΔt t
X , X vektörünün tt,t ve tt zamanlarındaki değerleri; t ise küçük ve sabit bir zaman adımıdır. t
X ve t
X in yeni durumları t anı için aşağıda verilmiştir. ) ( t 2 1 Δt t Δt t X X Xt ) 2 ( ) Δt ( 1 Δt Δt 2 t Xt Xt Xt X (4.7) (4.8) (4.7) ve (4.8) deki denklemlerdeki hata ( t )2 ile orantılıdır. Buna göre, t yarıya düşürülürse hata dörtte birine düşer.
Doğrusal problemler için tt zamanındaki yer değiştirmeyi bulabilmek için t zamanındaki hareket denkleminin ayrıklaştırılması gerekmektedir.
t t t t CX KX R X M (4.9) t X ve t
X (4.7) ve (4.8) denklemleri, (4.9) ile verilen hareket denkleminde yerine konursa, tamamiyle zaman içinde ayrıklaştırılmış bir hareket denklemi elde edilir.
2 t Δt t 2 t 2 t Δt Δt 2 1 Δt 1 Δt 2 Δt 2 1 Δt 1 X C M X M K R X C M (4.10)
Yeni oluşan denklem (4.10)‟da gösterilmiştir ve bu denklemden Xtt kolaylıkla
çözülebilir.
4.3.2 Merkezi Farklar Yöntemi Sonucunda Ortaya Çıkan Denklemin Çözümü
Δt t
X in çözümü için (4.10) nolu denklemde X , t XtΔt yer değiştirmelerine ihtiyaç vardır. X0 ve 0
X ın bilindiğini kabul edelim (X0ve 0
X biliniyorsa bunlar t=0 anı için (4.9) nolu denklemin içine yazılarak 0
X elde edilebilir). XtΔt yer değiştirmesini bulmak için ise (4.7) ve (4.8) nolu denklemleri kullanabiliriz. Bu denklemlerden birinde XtΔt yalnız bırakılıp diğer denklemde yerine yazılırsa XtΔt bulunabilir. Δt t X = t t t X X X 2 Δt Δt 2 (4.11)
t=0 anı için denklem (4.11)
Δt t X = 0 0 0 X X X 2 Δt Δt 2 (4.12) olur.
Denklem (4.11)‟de XtΔt nin son hali görülmektedir. (4.12) denkleminde ise t=0 anı için denklem tekrar oluşturulmuştur. (4.11) ve (4.12) ile verilen XtΔt ve X t ifadeleri (4.10) denkleminde yerine konursa, XtΔt‟nin değeri elde edilir. Bir sonraki prosedürün başlangıç koşulu olarak XtΔt=X ve t X =t XtΔt alınır.
4.4 Kapalı Zaman-Entegrasyon Prosedürü
4.4.1 Newmark Yöntemi
Newmark yöntemi kapalı uygulamalarda çok yaygın olarak kullanılan bir yöntemdir. Newmark yöntemi aşağıda gösterilen formülasyonları temel alır.
Δt Δ Δt t t t t X X X X 2 Δt t Δ (Δt) 2 1 Δt Xt Xt Xt Xt X t t X X X t Δt Δ (4.13) (4.14) (4.15) Eğer (4.15) numaralı denklem (4.14) ve (4.13) nolu denklemler yerine konulursa (4.16) ve (4.17) nolu denklemler elde edilirler.
1
t Δt Δt Δ t X X X X t t t 2 Δt t Δt t (Δt) 2 1 Δt X X X X X t t t (4.16) (4.17)Bu denklemlerde bulunan ve katsayıları çözümün kararlılığını sağlar. (4.17) no‟lu denklemden Δt t X çekilirse tΔt X = Xt Xt Xt X 1 2 1 Δt 1 ) Δt ( 1 ) Δt ( 1 2 Δt t 2 (4.18)
(4.18) no‟lu denklem elde edilir. Bu denklem (4.16)‟da yerine konursa
t t X X X X X 1 Δt 2 1 Δt Δt t Δt Δt t t (4.19)
denklemi bulunur. tt zaman adımındaki hareket denklemi
t Δt tΔt tΔt tΔt CX KX R X M (4.20)
Genel hareket denkleminde bulunan tt X , tt X ivme (4.18) ve hız (4.19) bilinmekteler. Bunlar (4.20)‟de yerine yazılır.
t t t t t t t t t X X X M X X X C R X M C K t 1 2 1 Δt 1 ) Δt ( 1 1 2 Δt 1 Δt ) Δt ( 1 Δt 2 2 4.21)
denklemi elde edilir. Newmark tarafından lineer ivmelenmeli durumlar için =1/4 ve =1/2 olarak önerilmiştir.
Başlangıç koşulu olarak da X ve 0 0
X yer değiştirme ve hıza ihtiyaç vardır. t=0 anı için (4.9) nolu denklemden de 0
X ivmesi elde edilmektedir.
4.4.2 Newmark Yöntemi Sonucunda Ortaya Çıkan Denklemin Çözümü
Δt t
X in çözümü için (4.21) nolu denklemdeki X , t t
X ve t X başlangıç şartlarının bilinmesi gerekmektedir. X0 , 0
X ‟ı başlangıç koşulu olarak alalım. X0ve 0
X biliniyorsa, 0
X ‟ı bulmak için bu değerler t=0 anı için (4.9) no‟lu denklemde kullanılabilir. X0, 0
X ve 0
X denklem (4.21)‟de kullanılarak X0Δt elde edilir.
Δt 0 X ,X0 , 0 X ve 0
X bilinenleri (4.19) denklemine konularak 0Δt
X elde edilir. Δt 0 X bulmak için X0Δt,X0 , 0 X ve 0
X bilinenleri (4.18) denkleminin içine yazılırlar ve 0Δt
X elde edilir. Bir sonraki prosedürün başlangıç koşulları olarak t X =XtΔt, t X = tΔt X ve t X = tΔt X alınır.
5. SAYISAL ÖRNEKLER
Yazılan programın doğru çalışıp çalışmadığı bu bölümde kontrol edilmektedir. Bu bölümde değişik kaynaklardan alınan örnekler yazılmış olan programla karşılaştırılmaktadır. En son örnekte, programa ek olarak yoğunluk da fonksiyonel olarak eklenmiş ve özgün bir örnek yapılmıştır.
5.1 Çubuğun Rijit Duvara Çarpması
Bu örnek, programın oluşturulma aşamasında programın doğru çalışıp çalışmadığını kontrol etmek amacıyla incelenmiştir. Kullanılan örnek [4] nolu referanstan alınmıştır. Bu örnekte düzlemsel birim uzama halindeki bir dikdörtgen çubuğun hızla rijit duvara vurması sonucu oluşan zamana bağlı yer değiştirmeler çözülmüştür. Referans [4]‟de çözüm yöntemi olarak Newmark seçilmiştir. Bu sayede bizim kullanmış olduğumuz Newmark yöntemini yazarın çözümüyle karşılaştırma ve Merkezi Farklar yöntemiyle aralarındaki farkları görme şansımız olacaktır. Şekil 5.1‟de çözülmüş olan problemin fiziksel halinin nasıl olduğuna dair bilgi vermektedir.
Şekil 5.1 Düzlem birim uzama halindeki bir çubuğun duvara çarpması. V
a
h V
Şekil 5.1‟de a*h boyutlarındaki bir çubuğun V hızıyla giderken önündeki rijit duvara vurması gösterilmektedir.
Bu örnek yazılmış olan sonlu elemanlar programıyla çözülürken, sonlu elemanlar modeli üzerine şu sınır şartları verilmiştir: Duvara çarpma anında çubuğun her noktasındaki hızı V‟dir. Çubuk duvara çarptıktan sonra çubuğun duvarla temas eden yüzeyinde hiçbir yer değiştirme olmayacaktır. Bu nedenle sonlu elemanlar modeli üzerindeki bütün düğüm noktalarına t=0 anı için V hızı verilmiş, çubuğun duvarla temas ettiği yüzeydeki düğüm noktaları için ise yer değiştirme “0” olarak alınmıştır. Şekil 5.2‟de çubuğun çözüm ağı gösterilmektedir. Toplam 21 elemandan oluşmakta ve 44 düğüm noktası bulunmaktadır. Yer değiştirmenin (duvara çarpma yüzeyindeki düğüm noktaları) sıfır olduğu düğüm noktaları 43 ve 44 üncü düğüm noktalarıdır ve şekil 5.2‟de gösterilmektedir. Şekil 5.5‟de gösterilen yer değiştirme zaman grafikleri 42‟nci düğüm noktası olarak alınmıştır. 42‟nci düğüm noktasının kullanılmasının nedeni, referans [4] deki yer değiştirme zaman grafiklerinin 42‟nci düğüm noktasına göre çizdirilmesidir. Analizlerde düğüm noktası kullanılarak sonuçların referans [4] ile karşılaştırılması yapılmıştır. Modelleme ve analizlerde kullanılan çubuğa ait bazı özellikler tablo 5.1‟de verilmiştir.
Şekil 5.2 Çubuğun çözüm ağı ve düğüm noktaları numaraları.
Tablo 5.1 Modellemede kullanılan dikdörtgen çubuğun özellikleri.
a 10.5 mm
h 0.5 mm
Elastisite modülü 100 Mpa
Yoğunluk 0.01 ton/mm3
Poisson oranı 0
Çubuğun dalga hızı 100 mm/sn
V 1 mm/sn
Eleman sayısı 21
Gama 0.5
Beta 0.25
t 0.00353
Merkezi farklar yönteminde kullanılan değişken
t 0.001765
Newmark ve Merkezi yöntemler kullanılarak çözülen problemin, 42‟nci düğüm noktasının yer değiştirme değeri, referans alınan kitaptaki değerlerle karşılaştırmalı olarak Tablo 5.2‟de verilmiştir.
Tablo 5.2 Bulunan yer değiştirme değerleri ve referans [4] den alınan değerlerin farkları.
Saniye Referans [4] Newmark Fark%
Merkezi
Farklar Fark % 3,53E-03 3,07E-03 3,01E-03 1,95 3,14E-03 2.28 7,06E-03 4,78E-03 4,59E-03 3,97 4,37E-03 8.57 1,06E-02 5,07E-03 4,82E-03 4,93 4,66E-03 8.08 1,41E-02 4,88E-03 4,72E-03 3,28 4,76E-03 2.459 1,77E-02 4,91E-03 4,76E-03 3,05 4,76E-03 3.05 2,12E-02 5,06E-03 4,80E-03 5,14 4,77E-03 5.731 2,47E-02 5,04E-03 4,76E-03 5,56 4,76E-03 5.555 2,82E-02 4,70E-03 4,76E-03 -1,28 4,77E-03 1.48 3,18E-02 4,96E-03 4,78E-03 3,63 4,77E-03 3.83 3,53E-02 5,02E-03 4,77E-03 4,98 4,77E-03 4.98
Şekil 5.4 42. düğüm noktasındaki yer değiştirmeler.
Şekil 5.4 çubuğun 42‟nci düğüm noktasının yer değiştirmesi, çarpma anından itibaren belirli bir saniyeye kadar Newmark ve Merkezi Farklar yöntemleri kullanılarak göstermektedir. Şekil 5.3 ise referans [5]‟den alınmıştır, 42‟nci düğüm noktasında meydana gelen yer değiştirmelerin grafiksel gösterimini vermektedir.
Şekil 5.5 Merkezi Farklar Yöntemi kullanılarak çubuk çarpmadan 1 saniye sonraki yer değiştirmeler.
Şekil 5.6 Newmark Yöntemiyle hesaplanan, çubuk çarpmasından 1 saniye sonraki yer değiştirmeler.
Şekil 5.5, Şekil 5.6 çubuğun çarpmadan 1sn sonra çubuk üzerindeki yer değiştirmeleri göstermektedir. Şekil 5.5‟te gösterilen değerler Newmark yöntemiyle bulunmuştur. Şekil 5.6‟da kullanılan yöntem Merkezi Farklar yöntemidir.
5.2 Silindirik Çubuk Üzerine Gelen Anlık Bir Kuvvet Sonucu Oluşan Yer Değiştirmeler
Burada yapılmış örnek ise Tanrıöver ve Şenocak‟ın [6] makalesinden alınmıştır. Makalede çözümler Newmark ve Merkezi farklar yöntemleri kullanılarak yapılmıştır. Bu makalenin seçilme nedeni, deneysel veriler içermesi ve her iki yönteminde kullanılmasındandır. Bu sayede yapılmış programın çözümleri ile makale çözümleri karşılaştırılabilmiş ve sonuçların doğruluğu kontrol edilmiştir.
Şekil 5.7‟de çözülen problemin fiziki hali hakkında fikir vermektedir. h*r boyutlarındaki silindirik çubuk alt kısmından sadece z yönünde mesnetlenmiştir. t=0 anında basamak fonksiyonu şeklinde silindirin üst yüzeyine basınç uygulanmıştır (F=2068427*H(t)). Buradaki H(t) Heaviside basamak fonksiyonudur.
Şekil 5.7 Silindirik bir çubuğun üzerine gelen anlık yükler. h
Şekil 5.8 Çubuğun yarısının çözüm ağı ve incelenen düğüm noktasının yeri.
Şekil 5.8‟de çubuğun yarısının çözüm ağı görülmektedir. Ayrıca aynı şekilde ileride incelenecek olan düğüm noktasının yeri gösterilmektedir. Bu düğüm noktasının seçilmiş olmasının sebebi, makalede yapılan çalışmalarda bu düğüm noktasının üzerindeki yer değiştirmelerin verilmiş olmasındandır. İlk olarak çubuk 5*100 elemandan oluşacak şekilde çözüm ağı oluşturulmuştur, eleman sayısı ve zaman adımı değiştirilerek sonuçların belirli bir değere yakınsaması sağlanmıştır. Analiz ve modelleme ile ilgili ek bilgiler Tablo 3‟de verilmiştir.
Tablo 5.3 Modellemede ve analizde kullanılan silindirik çubuğa ait özellikler.
Çubuk malzemesi 24S-T Al
r 12.7mm
h 254mm
Elastisite modülü 72395 Mpa
Yoğunluğu 2.7135 E-9 ton/mm3
Poisson oranı 0.3125
Çubuğun dalga hızı 5.1654 mm/ms
t=0 anında çubuğun üstüne uygulanan anlık basınç
2068427 H(t) Pa Newmak yönteminde kullanılan değişkenler
Gama 0.5
Beta 0.25
Ek A‟da çıkan sonuçlar gösterilmiştir. Şekil A.1, Şekil A.2, Şekil A.3 ve Şekil A.4‟de Newmark ve Merkezi Farklar kullanılarak, 0.008 milisaniye sonunda yapı üzerinde meydana gelen eksenel ve radyal yer değişimleri gösterilmektedir.
Şekil A.5 referans [6]‟dan alınmıştır ve deney sonucunda yapı üzerinde z=25.4 mm ve z=50.8mm deki düğüm noktalarına ait yer değiştirmeleri göstermektedir. Şekil A.6 ve Şekil A.7 sırasıyla Merkezi Farklar ve Newmark çözümleridir. Her bir grafikte üç değişik zaman adımında ve z=25.4mm, z=50.8mm deki düğüm noktalarının yer değiştirmeleri verilmiştir. Çıkan sonuçlar kendi içlerinde ve deney sonuçlarıyla uyumludur.
Newmark ve Merkezi Farklar yöntemlerini hız açısından karşılaştırmak için, analiz yapılırken, işlemci (CPU) zamanları kaydedilip Tablo 5.4‟de gösterilmiştir. Newmark yönteminin Merkezi Farklar yönteminden daha uzun sürmesinin nedeni, Newmark‟da yapılan işlem sayısının daha fazla olmasından kaynaklanmaktadır. Kullanılan bilgisayar Celeron 2.7 işlemciye ve 128 MB RAM‟e sahiptir
Tablo 5.4 Analiz zamanları (Cpu zamanı). ELEMAN
SAYISI
MERKEZİ FARKLAR NEWMARK
t
(sn) t(sn)
1E-7 8E-8 6E-8 1E-7 8E-8 6E-8
100X5 12.703125 26.0625 48.25 15.828125 34.625 63.7656 120X8 85.953125 160.453128 283.484375 110.078125 165.34375 411.187 150X10 90.90625 159.421875 297.734375 114.53125 210.375 422.828
5.3 Elastisitesi Değişen Silindirik Bir Plağın Statik Durumdayken Analitik Çözümleriyle Yazılan Program Sonuçlarının Karşılaştırılması
Burada yapılmış örnek ise J.N.Reddy‟nin makalesinden [2] alınmıştır. Bu kısma kadar elastisite sabit alınmıştı. Bu çalışmada ise statik halde üzerine “P” yayılı yükü, elastisitesi bir fonksiyona (5.1 denklemi) bağlı ve h kalınlığı boyunca değişen, kenarlarından mesnetlenmiş dairesel bir plağın statik analizi yapılacaktır ve çıkan sonuçlar N.Reddy „in makalesindeki [2] sonuçlarla karşılaştırılacaktır. Şekil 5.9‟da ise geometrinin özellikleri ve uygulanan basınç gösterilmiştir.
Şekil 5.9 Uygulanan kuvvet ve geometrinin şekli.
Modellemede kullanılan silindirik plağa ait özellikler.
n c n m h z h E h z h E z E 2 2 1 2 2 ) ( (5.1) z r
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 z E (z) /E c n=0 n=2 n=4 n=6
Şekil 5.10 Malzemenin kesit boyunca, elastisitesinin n ile değişimi.
(5.1) denkleminde elastisite modülünün değişim fonksiyonu verilmiştir. Elastisite özellikleri ise Em=0.8316 Gpa, Ec=2.1 Gpa‟dır Reddy, Wang ve Kitipornchai [2].
Yapının üst tarafı metal, alt tarafı ise seramiktir. Şekil 5.10‟da “n” sayısının değişimi sonucu elastisitenin kalınlık boyunca nasıl değiştiği gösterilmiştir.
Yarıçap (a) bütün analizlerde 1m olarak alınmıştır. h/a=0.05, h/a=0.15, h/a=0.2 oranları için ve n in 2, 4, 6, 8, 10, 15, 20, 25 değerleri alınarak analizler yapılmıştır. Yapılan analizlerde eleman sayısı 20x6 olarak seçilmiştir. Ancak h/a=0.2 durumunda 50*10 ve 100*20 eleman sayıları için de tekrar bir analiz yapılmış ve çözümün yakınsadığı gösterilmiştir.
h/a=0.05 için bulunan sonuçlar Tablo 5.5‟da gösterilmiştir. Tablo 5.5 h/a=0.05 için orta noktanın yer değiştirmesi.
Analitik çözüm(kabuk) Sonlu elemanlar çözümü Fark % n 1.67E-05 1.61E-05 3.94 0 9.18E-06 8.76E-06 4.61 2 8.39E-06 8.00E-06 4.73 4 7.99E-06 7.62E-06 4.64 6 7.73E-06 7.39E-06 4.47 8 7.56E-06 7.23E-06 4.43 10 7.29E-06 7.00E-06 4.07 15
7.15E-06 6.88E-06 3.81 20
7.05E-06 6.80E-06 3.61 25
h/a=0.15 için bulunan sonuçlar Tablo 5.6‟de gösterilmiştir. Tablo 5.6 h/a=0.15 için orta noktanın yer değiştirmesi.
Analitik çözüm(kabuk) Sonlu elemanlar çözümü Fark % n 6.74E-07 6.49E-07 3.8 0 3.67E-07 3.50E-07 4.77 2 3.35E-07 3.20E-07 4.69 4 3.19E-07 3.06E-07 4.30 6 3.10E-07 2.99E-07 3.56 8 3.03E-07 2.94E-07 3.05 10 2.93E-07 2.88E-07 1.72 15 2.87E-07 2.85E-07 0.78 20 2.83E-07 2.84E-07 0.1 25
h/a=0.2 için bulunan sonuçlar tablo 5.7‟de gösterilmiştir. Tablo 5.7 h/a=0.2 için orta noktanın yer değiştirmesi.
Analitik çözüm(kabuk) Sonlu elemanlar çözümü Fark % n 3.04E-07 2.96E-07 2.91 0 1.65E-07 1.59E-07 3.68 2 1.50E-07 1.45E-07 3.82 4 1.43E-07 1.39E-07 3.27 6 1.39E-07 1.35E-07 3.15 8 1.36E-07 1.34E-07 1.78 10 1.31E-07 1.31E-07 0.7 15 1.29E-07 1.30E-07 0.4 20 1.27E-07 1.29E-07 0.83 25
Şekil 5.11 Denklem (5.1)‟in değişik “n” lerde, eleman sayısına göre yakınsaması.
Plak boyutu 200x1000mm olarak alınmış, denklem (5.1)‟deki n‟nin farklı değerleri için, çözüm ağı oluşturulmuştur. Şekil 5.11‟de n‟nin değişik değerleri için hesaplanan birim uzamalar eleman sayısına bağlı olarak verilmiştir.
5.4 Elastisesi Değişen Dairesel Bir Plağın Üzerine Gelen Basınç Etkisiyle Plak Üzerindeki Bir Noktanın Yer değiştirmesinin Zamana Bağlı Bulunması
Şimdiye kadar yapılan bütün analizler statik olarak çözülmüştü. Bu örnek yazılmış olan programın dinamik analizle karşılaştırılmasını içermektedir. Araştırmalar sonucunda elastisesi değişen bir plağın dinamik analizlerine rastlanmadı. Bu nedenle analitik olarak bir örnek çözülmüş ve program sonuçlarıyla karşılaştırılmıştır.
z r w r 2 2 (5.2) r z r w (5.3)
r rr z E 2 1 ) ( (5.4)
r
z E 2 1 ) ( (5.5)Denklem (5.2), (5.3) silindirik koordinatlarda birim uzamaları, denklem (5.4), (5.5) silindirik koordinatlarda gerilmeleri göstermektedir. (5.2) ve (5.3) denklem (5.4) ve (5.5)‟de yerine yazılırsa.
z r r w r w z E rr 1 1 ) ( 2 2 2 z r w r r w z E 2 1 22 1 ) ( (5.6) (5.7)
(5.6) ve (5.7) denklemleri elde edilirler.
n c n m h z h E h z h E 2 2 1 2 2 E(z) (5.8)
Denklem (5.6) ve (5.7)‟de elastisite “z” ye bağlı olarak verilmiştir. Bu fonksiyon denklem (5.8)‟de görülmektedir Reddy, Wang ve Kitipornchai [2].
2 / 2 / h h rr r zdz M
/2 2 / h h zdz M (5.9) (5.10)Denklem (5.9) ve (5.10)‟da momentler görülmekte. Denklem (5.6),(5.7) ve (5.8) , denklem (5.9) ve (5.10)‟da yerine yazılırsa denklem (5.11) elde edilir.
2 / 2 / 2 2 2 2 1 1 1 2 2 1 2 2 h h n c n m r z dz r r w r w h z h E h z h E M (5.11)(5.1) denklemine bakıldığında Em =0.8316 Gpa Ec=2.1 Gpa n=2 v=0.288 alınmıştır. r r w r w D Mr 2 1 2 (5.12)
(5.11) denkleminin son hali (5.12)‟dir.
Aynı şeyler M içinde uygulanırsa denklem (5.13) elde edilir.
1 22 r w r r w D M (5.13)
D sabiti ise denklem (5.14)‟de gösterilmiştir.
2 3 1 0632 . 0 E h D c (5.14)
Referans [12] den silindirik bir plağın dinamik denge denklemleri alınmıştır.
2 2 4 t w P w D (5.15)
Malzemenin elastisitesinin kalınlık boyunca değişmesi (5.15) deki denklemin sadece “D” katsayısında değişikliğe sebep olur. Bundan dolayı da denklemin çözümünde bir değişiklik meydana gelmez. Bu denklemin çözümü referans [12]‟de ayrıntılı olarak verilmiştir.
Şekil 5.12 Analitik olarak çözülen geometri modeli.
Şekil 5.12‟de üzerine yük gelen plağın şekli verilmiştir. Bu geometriye ait özellikler tablo 5.9‟da verilmiştir. Yapı kenarlarından mesnetlenmiştir.
Tablo 5.8 Plağın özellikleri.
a 300mm b 48mm h(kalınlık) 15mm Ec 2.1e5 N/mm2 D 4.483 e7 p 0.000001 N/mm2 8e-9 ton/mm3
Şekil 5.14 Yapının merkezindeki yer değiştirmeleri
Şekil 5.15 Newmark ve Merkezi Fark yöntemi kullanılması sonucu yapının merkezinde oluşan yer değiştirmeler.
Yapının çözüm ağı ve incelenen düğüm noktası şekil 5.13‟de gösterilmiştir. Şekil 5.14‟de ise analitik çözümler, Şekil 5.15‟de ise Newmark ve Merkezi Fark sonuçları verilmiştir. Sonlu eleman sonuçlarının ve analitik çözümden çıkan sonuçların birbirleriyle örtüştükleri görülmektedir.
5.5 Elastisitesi Değişen Silindirik Bir Plağın Üzerine Gelen Anlık Bir Basınç Durumundaki Halinin İncelenmesi
Bu analizde diğer analizlerden farklı olarak, yoğunluk değişimi de hesaba katıldı. Yoğunluk, elastisitenin tanımlanmış olduğu gibi bir fonksiyonla ifade edilir.
Modellemede kullanılan silindirik plağa ait olan özellikler, aşağıda gösterilmiştir.
n h z b E z E 1 ) ( (5.15) n h z a z 1 ) ( (5.15)
En alttaki malzeme seramiktir. Es=420 (GP) s = 3200 kg/m3. Fonksiyonlardaki özellikler; Elastisite için n=1 b=1 yoğunluk için n=0.608 a=1
plağın en üst kısmında oluşan değerler E=210(GP) ve = 2100 kg/m3‟dir. Eksen boyunca malzeme özelliklerinin değişimi ise Şekil 5.16‟da verilmiştir.
0 0.5 1 1.5 2 2.5 -15 -10 -5 0 5 10 15 Elastisite E(z)/E yogunluk q(z)/q
Şekil 5.16 Plağın eksen boyunca malzeme değişimi.
Plağın geometrik özellikleri ve uygulanan kuvvet Şekil 5.17‟de gösterilmiştir.,
Şekil 5.17 Plak geometrisi.
Şekil 5.17 de gösterilen plağa ait fiziksel özellikler aşağıdaki gibidir. d=200mm
h=24,5mm a=6mm
uygulanan kuvvet P=3 (GPa)‟dır. Ek A‟da çıkan sonuçlar verilmiştir. Şekil A.8, Şekil A.9, Şekil A.10 ve Şekil A.11‟de Newmark ve Merkezi Farklar kullanılarak, 1e-4 saniye sonunda yapı üzerinde meydana gelen eksenel ve radyal yer değiştirmeleri gösterilmektedir. Ayrıca r=10 mm ve z=254 mm deki noktanın eksenel ve radyal yer değiştirmeleri için değişik zaman adımlarında analizler yapılmıştır ve Şekil A.12, Şekil A.13, Şekil A.14, Şekil A.15‟de verilmiştir.
6.SONUÇLAR VE TARTIŞMA
Bu çalışmada özellikleri fonksiyonel olarak değişen malzemelerin özellikleri ve mekanik davranışları incelendi. Bu malzemelerin dinamik yükler altındaki deformasyon alanlarının belirlenmesi için gerekli sonlu elemanlar formülasyonları elde edildi. Açık ve kapalı yöntemler kullanılarak elde edilmiş formülasyonlar çeşitli problemlere uyarlanarak sonuçlar karşılaştırıldı. Bu çalışma boyunca nispeten basit sınır koşullarına sahip ve özellikleri fonksiyonel olarak değişen yapıların dinamik davranışları incelenmiştir. Yapının mukavemet özelliklerinin zamanla değişmediği ve sistemin tamamen doğrusal olduğu farz edilmiştir.
Silindirik bir çubuk ve plak için yapılan sonlu elemanlar analizleri göstermektedir ki Newmark Yöntemi, Merkezi Farklar Yöntemine göre daha iyi sonuçlar vermiştir. Belirli bir hata sınırı için Newmark Yöntemi daha geniş zaman adımlarının kullanılmasını ve böylece sayısal çözümün çok daha kısa sürede sonuçlanmasını sağlamaktadır. Ancak Tablo 5.4 incelendiğinde, aynı zaman adımı kullanılarak yapılan analizlerde merkezi farklar yönteminin daha hızlı sonuç verdiği görülmektedir. Bu nedenle Newmark ve benzeri kapalı algoritmalar büyük zaman adımlarının kullanılabileceği ve yüksek hassasiyet gerektiren analizler için elverişlidir. Merkezi Farklar ve benzeri açık yöntemler ise çarpışma, patlama gibi çok ani gelişen ve çok küçük zaman adımlarının kullanılmasını gerektiren fiziksel olayların modellenmesi için elverişlidir.
Bütün bölümlerde doğrusal sistemlerin modellenmesi için yapılan formülasyonlar doğrusal olmayan modellere genişletilebilir. Malzeme, geometri, yük ve kontak modeliyle oluşan doğrusal olmayan davranışın katılık, yük ve sönüm matrisine olan etkileri incelenebilir. Daha farklı açık ve kapalı yöntemlere ait algoritmaların performansları benzeri problemler için karşılaştırılabilir. Günümüzde artık ön plana çıkan belirli bir fiziksel sistemin modellenmesi için maliyeti, tasarım zamanı ve çözüm doğruluğu optimize edilmelidir. Bu nedenle pahalı olan deneysel yöntemler yerine açık ve kapalı yöntemlerin birlikte kullanıldığı hibrit yöntemler geliştirilebilir.
KAYNAKLAR
[1] Chiu, T. -C. and Erdogan F., 1999. One-Dimensional Wave Propagation in A Functionally Graded Elastic Medium, Journal of Sound and Vibration, 222(3), 3-487.
[2] Reddy, J. N., Wang, C. M. and S. Kitipornchai, S., 1999. Axisymmetric
Bending of Functionally Graded Circular And Annular Plates, Eur, J, Mech. A/Solids, 18, 185-199.
[3] Dokainish, M. A. and Subbaraj, K., 1989. A Survey of Direct Time-Integration Methods in Computational Structural Dynamics-1. Explicit Methods, Computers Structures, 32, 1371-1386.
[4] Dokainish, M. A. and Subbaraj, K., 1989.A Survey of Direct Time-Integration Methods in Computational Structural Dynamics-2. İmplicit Methods, Computers Structures, 32, 1387-1401.
[5] Griffths, D. V. and Smith, I. M., 1997. Programing the Finite Element Method,
Chichester: Wiley.
[6] Tanrıöver, H. ve Şenocak, E., Çubukta Gerilme Dalgalarının İlerlemesi Problemleri için Zaman İntegrasyonu Yöntemlerinin Karşılaştırılması, TUMTMK XIII ulusal Mekanik kongresi, Gaziantep 8-12 Eylül-2003.
[7] Li, Y., Ramesh, K.T. and Chin, E. S. C., 2001. Dynamic Characterization of Layerd and Graded Structures Under İmpulsive Loading, İnternational Journal of Solids and Structures, 38, 6045-6061. [8] Berezovski, A., Engelbrecht, J. And Maugin, G. A., 2003. Numerical
Simulation of Two-Dimensional Wave Propagation in Functionally Graded Materials, European Journal of Mechanics A/Solids,22,257-265.
[9] Gazonas, G.A.,Thamburaj, P. and Santare, M. H., 2003.The Use of Graded Finite Elements in the Study of Elastic Wave Propagation in Continuously Nonhomogeneous Materials, İnternational Journal of Solids and Strucures, 40, 5621-5624.
[10] Liu, G. R. and Han, X., 2002. Effects of SH Waves in a Functionally Grade Plate , Mechanical Reserch Communications, 29, 327-338.
[11] Velasco, V. R., Zarate, J. E., 2002. Elastic Wave in Graded-Composition System, Superlattices and Microstructures, 28, No:3.
[12] Meirovitch, Leonard., 1967. Analytical Methods in Vibrations, The Machillan Company Collier-Macmillan Limited, London.
EK A
Şekil A.1 Newmark yöntemi ile t=0.008 milisaniyedeki çubuğun üzerindeki eksenel yer değiştirmeler.
Şekil A.2 Merkezi farklar yöntemi ile t=0.008 milisaniyedeki çubuğun üzerindeki eksenel yer değiştirmeler.
Şekil A.3 Newmark metodu ile t=0.008 milisaniyedeki çubuğun üzerindeki radyal yer değiştirmeler.
Şekil A.4 Merkezi farklar metodu ile t=0.008 milisaniyedeki çubuğun üzerindeki radyal yer değiştirmeler.
Şekil A.5 z=25.4 mm ve z=50.8 mm deki yer değiştirmeleri göstermektedir [6].
Şekil A.6 Merkezi Farklar Yöntemi kullanılarak değişik zaman adımlarında ve değişik noktalardaki radyal yer değiştirmeler.
Şekil A.7 Newmark kullanılarak değişik zaman adımlarında ve değişik noktalardaki radyal yer değiştirmeler.
Şekil A.8 Newmark Yöntemi kullanılarak, 1e-4 zaman aralığında, 50x10 eleman sayısında, r=10 mm ve z=254 mm noktasındaki düğümün eksenel yer değiştirmesi.
Şekil A.9 Merkezi Farklar Yöntemi kullanılarak, 1e-4 zaman aralığında, 50x10 eleman sayısında, r=10 mm ve z=254 mm noktasındaki düğümün radyal yer değiştirmesi.
Şekil A.10 Merkezi Farklar Yöntemi kullanılarak, 1e-4 zamanında, 50x10 eleman sayısında, r=10 mm ve z=254 mm noktasındaki düğümün eksenel yer değiştirmesi.
Şekil A.11 Merkezi Farklar Yöntemi kullanılarak, 1e-4 zamanında, 50x10 eleman sayısında, r=10 mm ve z=254 mm noktasındaki düğümün eksenel yer değiştirmesi.
Şekil A.12 r=10 mm ve z=24.5 mm deki düğüm noktasının, Newmark Farklar Yöntemi kullanılarak 0.1 milisaniye zaman aralığındaki radyal yer değiştirmesi.
Şekil A.13 r=10 mm ve z=24.5 mm deki düğüm noktasının, Newmark Farklar Yöntemi kullanılarak 0.1 milisaniye zaman aralığındaki eksenel yer değiştirmesi.
Şekil A.14 r=10 mm ve z=24.5 mm deki düğüm noktasının, Merkezi Farklar Yöntemi kullanılarak 0.1 milisaniye zaman aralığındaki radyal yer değiştirmesi.
Şekil A.15 r=10 mm ve z=24.5 mm deki düğüm noktasının, Merkezi Farklar Yöntemi kullanılarak 0.1 milisaniye zaman aralığındaki eksenel yer değiştirmesi.
ÖZGEÇMİŞ
1979 yılında Ankara‟nın Güdül İlçesinde doğdu. İlk, orta ve liseyi Ankara‟da tamamladı.Liseye 1993 yılında İncirli Süper Lisesinde başladı ve 1997 yılında İstanbul Teknik Üniversitesi Uçak Mühendisliğini Kazandı. 2002 yılında mezun oldu. Bitirme Tezi olarak insansız hava araçlarının gövde tasarımı ve analizi aldı. 2002 yılında İstanbul Teknik Üniversitesi yüksek lisans bölümü Katı Cisimlerin Mekaniğine kabul edildi.