HibSis BİLGİSAYAR YAZILIMI - ANKARA ÜNİVERSİTESİ FEN BİLİMLERİ ENSTİTÜSÜ DOKTORA TEZİ MELEZ ( H

108

EK 1 ELASTİSİTE KURAMI VE ELASTİK DALGA DENKLEMİ

Fizikte dalga genel olarak iki gruba ayrılır: ilerleyen ve durağan dalgalar. Sismik dalgalarda bu iki dalga türündendir. İlerleyen dalgalar, sismik kaynaktan itibaren her doğrultuda ilerleyen dalgalardır. Durağan dalgalara aynı zamanda yerin serbest salınımları’da denilir ve yerin bir bütün olarak salınımını gösterir. Durağan dalgalar ancak büyük depremlerle üretilebilir.

Enerjinin mekansal odaklanması açısından bakıldığında, sismik dalgalar hacim (cisim) ve yüzey dalgaları olmak üzere iki farklı gruba ayrılır. Hacim dalgaları bulunduğu hacim içerisinde yayılırken, yüzey dalgaları ortam yüzeyi boyunca yayılma gösterirler.

Havada ilerleyen ses (akustik) dalgaları veya boşlukta ilerleyen elektromanyetik dalgalarda hacimsel dalgalara birer örnek oluşturur. Yüzey dalgalarına örnek ise suya atılan taşın oluşturduğu dalga gösterilebilir. Fakat su dalgaları ile sismik dalgalar arasındaki en temel fark, dalgayı oluşturan kuvvetlerdir. Su dalgalarını oluşturan en temel kuvvet, yerçekimi kuvvetidir (örneğin, yerçekimi kuvvetinin üst üste binmesi ve yerkürenin kendi ekseni etrafında dönmesinden ötürü oluşan merkezkaç kuvveti). Bu nedenle, su dalgalarına gravitasyonel dalgalar da denir. Buna karşın sismik dalgaları oluşturan kuvvetler ise ortamdaki elastik kuvvetlerdir. Bu noktada, gravite kuvvetleri ihmal edilebilecek kadar küçüktür ve çoğu durumda ihmal edilir. Sonuç olarak, sismik dalga yayılımının temeli elastisite teorisine dayanır. Bu kısımda, temel elastisite teorisi içerisinde gerilme (stress) ve yamulma (strain) ilişkisi ile sayısal çözüm yöntemlerinde kullanılan elastik dalga denklemi ifadeleri verilmiştir. Türetimlerde iki temel yaklaşımdan (Langrange ve Euler) Langrange türü yaklaşım benimsenmiştir.

Gerilme, uygulanan dış yükler altında, bir cismin her bir birim alanındaki içsel kuvvetlerin dağılımının bir ölçüsüdür. Gerilme analizinde amaç, verilen kuvvet ve koşullar altında elemanda oluşan yer değiştirmelerin veya deformasyonların belirlenmesidir. Eksensel gerilme,

= F σ

109

ile verilir (Aki ve Richards 1980, Lay ve Wallace 1995 ). Eğer kuvvet alana dik ise bu gerilmeye normal gerilme veya basınç gerilmesi, kuvvet alanın bir parçasına teğet ise kayma gerilmesi adını alır.

⎟⎟

⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

zz zy zx

yz yy yx

xz xy xx

σ τ τ

τ σ τ

τ τ σ

Bu bağıntı ile verilen dizeye “gerilme tensörü” denir (Terzaghi 1962). Tensör matematiksel olarak genelleştirilmiş bir değer veya tüm skaler, yöney, dizey ve doğrusal işleçleri içeren geometrik durumun özel bir halidir. Burada σ normal, τ makaslama gerilmelerini göstermektedir. Gerilme tensörü bir küpün altı yüzüne uygulanan gerilmelerden oluşur. Gerilme için kullanılan ilk indis gerilmenin doğrultusunu, ikinci indis ise gerilmenin etki ettiği yüzeyi gösterir. Şekil 1’de gerilme bileşenlerinin etkili olduğu küp eleman görülmektedir.

Şekil 1 Küp eleman yüzeyindeki gerilme bileşenleri

xx dx

x σ + ^∂σ

∂

xy dx

x τ + ^∂τ

∂

xz dx

x τ + ^∂τ

yz ∂

yz dy

y τ + ^∂τ

∂

yy dy

y σ + ^∂σ

∂

zz zz dz

z σ + ^∂σ

∂

zy dz

z τ + ^∂τ

∂

zx dz

z τ + ^∂τ

∂

yx dy

y τ + ^∂τ

∂

110

Elastik bir cisim gerilme altında hacim ve şekil değişikliğine uğrar. Buna birim deformasyon ya da yamulma adı verilir. Yamulma, bir cismin gerilme uygulandıktan sonra başlangıca oranla şekil ve boyutlarında meydana gelen değişikliğin miktarı olarak tanımlanır. Bu biçim değişikliği, yamulma tensörü ile ifade edilebilir. Yamulma göreceli bir ölçümdür ve bu nedenle boyutsuzdur. Şekil 2’de l0 uzunluğundaki bir elemana eksensel gerilme uygulanırsa, eleman boyundaki değişim Δl= l−l₀ olur.

Elemanda oluşan eksensel deformasyon ise,

Δl ε=

ile verilir (Terzaghi 1962).

Şekil 2 Eksensel deformasyon

⎟⎟

⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

zz zy zx

yz yy yx

xz xy xx

ε ε ε

Yamulma ile ilişkili diğer bir nicelik ise yer değiştirmelerdir. Bunlar kısaca malzemenin herhangi bir noktası tarafından hareket ettirilen uzaklıktır. Malzemede bir yamulma oluştuğu zaman, yer değiştirme miktarları noktadan noktaya göre değişir.

111

Elastisite teorisi temelde kuvvetlere maruz kalan bir maddenin davranışının (deformasyonu veya büyük gerilmeler altında nihayi bozulmaların) incelenmesidir.

Günümüzün modern elastisite teorisinin temelleri Navier (1821) tarafından denge denklemlerinin ve elastik katıların hareketlerini göstermesiyle atılmıştır. Gerilme ve yamulma arasındaki ilişki madde türüne ve çok sayıda diğer etkenlere bağlıdır. Gazlar ve sıvılar için bu ilişki katılardan oldukça farklıdır. Hatta aynı fazdaki katı maddeler arasında dahi belirgin gerilme-yamulma farklılıkları görülür. Maddenin gerilme yamulma ilişkisini reoloji bilimi inceler. Katı maddelerdeki gerilme-yamulma ilişkisinin oldukça karmaşık olmasından dolayı genel olarak yapılan modellerde basitleştirme işlemine gidilir. Maddeler reolojik davranışlarına göre temel üç grupta incelenir. Bunlar;

elastik, viskoz ve plastik davranışlarıdır. Üzerinde kuvvet uygulanan maddelerde kuvvetin kaldırılmasıyla kalıcı yamulmanın oluşmadığı maddelere elastik maddeler denir. Bu tür davranış gösteren maddelerde gerilme-yamulma arası ilişki doğrusaldır.

σ =cε

Madde üzerine küçük bir yük uygulandığında zamanla artan bir yamulma oluşursa bu tür maddelerde gerilme-yamulma ilişkisi zamana bağlıdır ve Newton türü maddeler denir. Bu tür maddelerde gerilme, yamulma oranı ile doğrusal ilişki gösterir

c t σ = ^∂ε

∂

ve viskoz türü davranışı temsil eder. Madde üzerine uygulanan kuvvet belirli bir sınır değeri aştığında madde kuvvetin kaldırılmasından sonra da sıvı gibi akmaya başlar. Bu davranışa plastik davranış denir. Yüksek basınç altındaki birçok metal ve makro moleküllü maddeler plastik davranış sergilerler. Şekil 3’te sözü edilen üç tür davranış alanları gerilme-yamulma eğrisi üzerinde gösterilmiştir. Denge denklemlerinin çözülebilmesi için üzerine kuvvetlerin etkidiği bir hacim elemanında gelişen gerilme ile yamulma arası ilişkinin belirlenmesi gerekir. Elastik davranış bölgesi için geçerli olan Hooke ilkesi, yalnızca kuvvetlerin etkidiği doğrultularda oluşan yamulmayı tanımlar.

112

Bununla birlikte, gerilme ve yamulma tensör karekterinde ilişkili büyüklüklerdir. Bu nedenle, gerilme ve yamulma tensörleri arasında da doğrusal ilişki varsayımı yapılarak genelleştirilir (Love 1944, Terzaghi 1962).

Şekil 3 Katı maddelerde elastik, viskoz ve plastik davranış

1,3 1,3

ij ijk l k l ijk l k l

k l

c c

σ ε ε

= =

= ≡ ∑ ∑

(1.1)

(1.1) ifadesi genelleştirilmiş Hooke ilkesidir. Burada σ gerilmeyi, ε yamulmayı, c_ijkl ise elastik katsayı dizeyini gösterir. (1.1) ifadesi birçok maddenin, örneğin kristal yapıları, izotrop maddeleri ve birçok diğer anizotrop maddelerin davranışını açıklamakta oldukça başarılıdır. c elastik katsayı dizeyi genel durumda 3_ijkl ⁴=81 elemanlıdır. Gerilme ve yamulma tensörlerindeki bakışımlılık özelliğinden dolayı bağımsız katsayı adedi 6x6=36 olur. Elastik davranış gösteren maddeler için birinci ve ikinci indisler arası değişim özelliğinden (c_ijkl =c_klij) dolayı bağımsız katsayı adedi 21 olur. 21 adet elastik katsayı daha ziyade triklinik kristalografik yapılarda görülür.

Yüksek simetri özelliği gösteren yapılarda bağımsız katsayı adedi daha da azalır.

Örneğin, monoklinik yapılarda 13, rombik yapılarda 9 ve kübik yapılarda ise 3 bağımsız katsayı ile ifade edilir. İzotropik yapılarda tüm doğrultularda fiziksel özelliklerin aynı olmasından dolayı 2 bağımsız katsayı ile gösterilir. Bu iki katsayı λ ve μ ile simgelenir. Sonuç olarak, izotropik ortam için genelleştirilmiş Hooke ilkesiyle gerilme

113

yamulma ilişkisi (Novotny 1999, Lay ve Wallace 1995) Einstein toplama ilkesine uygun olarak:

( )

2 2

ij ijkl kl

ij kl ik jl il jk kl

kk ij ij ji

kk ij ij

ij ij

τ

ε

λδ δ μ δ δ δ δ ε λε δ μ ε ε

λε δ με λθδ με

⎡ ⎤

=⎣ + + ⎦

= + +

= +

(1.2)

(1.2) ifadesinde θ ε= ₁₁+ε₂₂+ε₃₃ hacimsel değişimi (dilatasyon) gösterir. δ _ij Kronecker deltasıdır (i=j için δ

ij =1, i ≠ j için δ

ij = 0 olması dönme hareketinin olmadığını yalnızca öteleme hareketinin olduğunu gösterir). λ ve μ katsayılarına Lamme katsayıları adı verilir. μ katsayısı makaslama modülünü gösterirken, λ katsayısının açık fiziksel bir anlamı yoktur.

Bir madde üzerine etkiyen toplam kuvvet ve toplam döndürme (tork) kuvveti sıfır ise madde statik dengededir. Sürekli ortamın bir parçası ele alındığında yamulmanın oluştuğu durumda da aynı koşul geçerlidir. Sözü edilen bu denge koşulunun (döndürme kuvveti hariç) integral ifadesi:

) 0

( +

∫∫∫

∫∫

S V

FdV dS

T ^υ (1.3)

yapısında yazılabilir. (1.3) ifadesinde V gözönüne alınan sürekli ortam parçasının hacmini, S, V hacmini saran yüzeyi, T^(ν) S yüzeyi normali (ν) doğrultusunda uygulanan kuvveti, F ise hacimsel kuvveti gösterir. Herhangi bir t anında sonsuz küçük hacmi V

=V(t), hızını v ve ρ’da yoğunluğunu gösterirse, sonsuz küçük V hacmindeki noktanın momentumu:

114 dV

v P

∫∫∫

= ρ (1.4)

yazılabilir. P momentum ifadesinin zamana bağlı türevi eylemsizlik (atalet) kuvvetini verecektir. (1.3) ifadesinin sağ tarafına (1.4) terimi eklenerek,

∫∫∫

∫∫

⁺ ⁼

V V

dt vdV FdV d

T⁽^{υ )} ρ (1.5)

elde edilir. (1.5) eşitliği, her bir t zaman değeri için madde üzerinde etkili olan hacim ve yüzey kuvvetleri toplamının, oluşan eylemsizlik kuvvetine eşit olması gerektiğini ifade eder ve sürekli ortam için hareket denkleminin integral ifadesidir. Elastik dalga yayılımı problemlerinde çoğu zaman sayısal hesaplamalardaki kolaylıklarından ötürü integral ifadesi yerine diferansiyel yapıdaki ifadeler tercih edilir. Bu nedenle, (1.5) ile verilen hareket denkleminin diferansiyel yapısını elde etmek için Gauss teoreminden yararlanılabilir:

∫∫∫

⁼

∫∫

V S

dS n F dV F

div G G

. )

( (1.6)

Burada FG

sürekli türevlenebilir kuvvet yöneyini, nG

S yüzeyi normalini gösterir birim yöneydir. x x x x x= ( , , )₁ ₂ ₃ zamandan bağımsız koordinatları ve div F( ) = df /dx _j _j özelliği kullanılırsa, (1.6) ile verilen Gauss ifadesi farklı yapıda yazılabilir:

∫∫∫

^∂ ⁼

∫∫

V S

j j j

j dV f n dS

f G G

. (1.7)

(1.3) ile verilen denge koşulu ifadesi i.terim için

115

) 0

( +

∫∫∫

∫∫

i S

i dS f dV

T^υ

(1.8)

fj=τ eşitliği ve (1.7) ifadesi kullanılırsa (1.8) ifadesindeki yüzey integrali hacim integrali yapısında yazılabilir:

∫∫∫

⎟⎟ ⁼

⎠

⎞

⎜⎜

⎝

⎛ +

∂

i j

ji f dV

x 0

τ (1.9)

(1.9) ifadesininin sıfıra eşit olması ancak integrantın sıfıra eşit olmasıyla mümkündür.

Diğer yandan integrantın dV hacmi içerisinde sürekli olduğu varsayılmasından dolayı:

∂ +

∂

i j

ji f

τ (1.10)

olmalıdır. (1.10) ifadesi denge koşulunun diferansiyel ifadesidir. Birim hacim elemanı için eylemsizlik kuvveti:

) t( v F_iner ρ

∂

− ∂

= (1.11)

ile verilir. Burada ρyoğunluk, v hız, t ise zamandır. Yoğunluğun zamanla değişimi ihmal edilebilir olduğu varsayılırsa, hızın (v) Lagrange koordinatlarındaki ifadesi

) , , ,

(x₁ x₂ x₃ t v

v= şeklindedir. x1, x2, x3 koordinatları ilksel pozisyonu ifade eder ve t zamanından bağımsızdır.

Zamana göre toplam türev:

116

2 2

t u dt

F_iner dv

∂

− ∂

−

= ρ ρ (1.12)

ile verilir. Burada u yerdeğiştirme yöneyidir. d’Alembert ilkesine göre, hareket denklemi denge koşulundan türetilebilir. (1.10) ile verilen denge koşulu ifadesine atalet kuvveti eklenirse:

2 2

i ij

i j

u f

t x

ρ^∂ = ^∂τ +

∂ ∂ (1.13)

ifadesi elde edilir. Bu ifade elastik dalga teorisinin temel denklemidir.

Birim hacim elemanının ivmesi, yer değiştirme vektörünün zamana göre ikinci türevi ile verilmektedir. Newton’un ikinci yasası yani hareket ilkeside (F = m.a) kullanıldığında (1.13) ifadesinde verilen hareket denklemi’ni yazmak mümkündür (Aki ve Richards 1980, Lay ve Wallace 1995, Moczo vd. 2004).

Eşitliğin sol tarafındaki terimler ortamdaki gerilme gradyentleri ve cisim kuvvetleriyle ilişkili yoğunluk ağırlıklı ivmeyi göstermektedir. Cisim kuvvetlerinin ( fi ) olmadığı durumda, tekdüze dalga denklemi elde edilir.

2 2

i ij j

t x

ρ^∂ = ^∂τ

∂ ∂ (1.14)

Bu eşitliği çözebilmek için, yer değiştirmeler (u) anlamında gerilmelerin (τ) bir ifadesi olarak gerilme-yamulma ilişkisinden (Hooke yasasından, denklem (1.1)) yararlanılır.

Kartezyen koordinatlarda bu hareket denklemi,

117 z

y t x

∂ +∂

∂

= ∂

∂

∂ ₁₁ ₁₂ ₁₃

2 1

2 σ σ σ

z y

x t

∂ +∂

∂

= ∂

∂

∂ ₂₁ ₂₂ ₂₃

2 2

2 σ σ σ

ρ (1.15)

z y

x t

∂ +∂

∂

=∂

∂

∂ ₃₁ ₃₂ ₃₃

2 3

2 σ σ σ

şeklinde olur (Virieux 1986). Burada u1,u2,u3 ifadeleri x,y,z yönlerindeki yerdeğiştirmelerdir. (1.2) ifadesinde verilen Hooke yasası, kartezyen koordinatlarda, yer değiştirme ile gerilmeler arasında 6 lineer denklemden oluşur. Hooke yasasından, gerilme ve yamulma arasındaki bağıntılar ise,

11 2 11

σ =λθ + με

22 2 22

σ =λθ+ με

33 2 33

σ =λθ+ με (1.16)

12 2με

σ =

13 2με

σ =

23 2με

σ =

şeklinde olurlar (Lay ve Wallace 1995). (1.15)’te verilen ifadelerdeki yamulma terimlerinin yerdeğiştirmeler türünden tanımlamaları ise,

x u

∂

= ∂ ¹ ε11

y u

∂

= ∂ ² ε22

z u

∂

= ∂ ³ ε33

1 2

i j ij

j i

u u

x x

ε = ^⎛⎜⎜⎝∂^∂ +^∂∂ ^⎞⎟⎟⎠

118

⎟⎟⎠

⎞

∂ +∂

⎜⎝

⎛

∂

= ∂

y u x u₂ ₁

12 2

ε 1 (1.17)

⎟⎟⎠

⎞

∂ +∂

⎜⎝

⎛

∂

= ∂

y u z u₂ ₃

23 2

ε 1

⎟⎠

⎞

∂ +∂

⎜⎝

⎛

∂

= ∂

z u x u₃ ₁

13 2

ε 1

şeklinde olur. (1.16) ve (1.17) denklemleri (1.15) denkleminde yerine yazılırsa 3-B kartezyen koordinatlarda hareket denklemi elde edilir.

⎜⎜

⎝

⎛

⎟⎟⎠

⎞

∂

∂ +∂

∂

∂ + ∂

⎟⎟+

⎠

⎞

∂

⎜⎜ ∂

⎝

⎛ +

∂ + ∂

∂ =

∂

z x

u y x

u z

u y

u x

u t

u ² ₂ ² ₃

2 1 2 2

1 2 2

1 2

) ( )

(λ μ μ λ μ

⎜⎜

⎝

⎛

⎟⎟⎠

⎞

∂

∂ +∂

∂

∂ + ∂

⎟⎟+

⎠

⎞

∂

⎜⎜ ∂

⎝

⎛ +

∂ + ∂

∂ =

∂

z y

u y x

u z

u x

u y

u t

u ² ₁ ² ₃

2 2 2 2

2 2 2

2 2

) ( )

(λ μ μ λ μ

ρ (1.18)

⎜⎜

⎝

⎛

⎟⎟⎠

⎞

∂

∂ +∂

∂

∂ + ∂

⎟⎟+

⎠

⎞

∂

⎜⎜ ∂

⎝

⎛ +

∂ + ∂

∂ =

∂

z y

u z x

u y

u x

u z

u t

u ² ₁ ² ₂

2 3 2 2

3 2 2

3 2

) ( )

(λ μ μ λ μ

Bu denklemlerde

ρ μ λ+2

p =

V ve

= μ

Vs şeklinde yazılıp düzenlenirse,

⎜⎜

⎝

⎛

⎟⎟⎠

⎞

∂

∂ +∂

∂

− ∂

⎟⎟+

⎠

⎞

∂

⎜⎜ ∂

⎝

⎛ +

∂ + ∂

∂

= ∂

∂

z x

u y x V u z V

u y

V u x V u t

s p s

p 3

2 2 2 2

2 2

1 2 2

1 2 2 2

1 2

) (

⎜⎜

⎝

⎛

⎟⎟⎠

⎞

∂

∂ +∂

∂

− ∂

⎟⎟+

⎠

⎞

∂

⎜⎜ ∂

⎝

⎛ +

∂ + ∂

∂

= ∂

∂

z y

u y x V u z V

u x

V u y V u t

s p s

p 2 3

1 2 2

2 2

2 2 2

2 2 2 2

2 2

)

( (1.19)

⎜⎜

⎝

⎛

⎟⎟⎠

⎞

∂

∂ +∂

∂

− ∂

⎟⎟+

⎠

⎞

∂

⎜⎜ ∂

⎝

⎛ +

∂ + ∂

∂

= ∂

∂

z y

u z x V u y V

u x

V u z V u t

s p s

p 2

2 1 2 2 2 2

3 2 2

3 2 2 2

3 2

) (

119

3-boyutlu Elastik Dalga Denklemi elde edilir. Bu denklemlerde V_p² =V_s² alınırsa,

⎟⎟⎠

⎞

∂ +∂

⎜⎜⎝

⎛

∂

= ∂

∂

2 1 2 2

1 2 2

1 2 2 2

1 2

z u y

u x

V u t

p (1.20)

şeklinde 3-boyutlu Akustik Dalga denklemi elde edilir.

120

EK 2 2B ELASTİK DALGA DENKLEMİNİN SÖNÜM DİZEYİ KULLANILARAK ÇÖZÜMÜ

2-Boyutta elastik dalga denkleminin SEY’ne göre çözümünün dizey ifadesi olan (4.28) denklemi, C sönüm dizeyi olmak üzere

M U CU+ + K U = F (2.1)

olarak ifade edilir. Bu dizey formunun Newmark yaklaşımı kullanılarak, t + tΔ zamanı için çözüm ifadeleri (α=0.5, β=0.25 olmak üzere) (Zienkiewicz ve Morgan 1983, Zienkiewicz ve Taylor 1991, Ulrich 2004)

1 2

a = 1

βΔt , a = ₂ t α

βΔ , ₃ 1 a =

βΔt , ₄ 1 a = - 1

β ^, ^{a =}⁵ ^{- 1} α β ^,

a = t - 2 2

α β

⎛ ⎞

Δ ⎜ ⎟

⎝ ⎠ , a = t 1 - 7 Δ

(

)

, a = ₈ αΔ t (2.2)

1 2 3 4 5 6 7 8

a , a ,a , a ,a ,a ,a ,a katsayılar olmak üzere, etkin (effective) sıkılık dizeyi, Kef,

ef 1 2

K = +K a M a C+ (2.3)

etkin kuvvet vektörü Fef

. .. . ..

ef 1 3 4 2 5 6

F (t+ Δ =t) f t( + Δ +t) M a U⎡⎢⎣ t+ a Ut+ a Ut⎤⎥⎦ +C a U⎡⎢⎣ t + a Ut+ a Ut⎤⎥⎦ (2.4)

K .efU_t_+Δ_t =F t_ef( + Δ t) (2.5)

[ ]

.. .

t+ t 1 t+ t 3 t 4 t

U Δ = a U Δ −U_t - a U - a U (2.6)

121

t t t+ t

t+ t = U + a U7 a U8

U Δ + Δ (2.7)

olarak ifade edilirler. (2.1) ifadesinin Merkezi Farklar yaklaşımı kullanılarak elde edilen çözüm ifadeleri için,

1 1

t 2 t t

t U U U

U t

+ − + −

≈ Δ

(2.8)

(2.9)

tanımlanan birinci ve ikinci türev eşitlikleri (2.1) ifadeisnde yerine yazılarak,

1 2 1 2 1.( . )

k k k k

U ⁺ = U −U ⁻ + Δt M⁻ Fs−K U (2.10)

(2.1) denkleminin SFY çözüm ifadesi elde edilir.

1 1

t t

t U U

U t

+ − −

≈ Δ

122

EK 3 GLOBAL-LOKAL KOORDİNAT DÖNÜŞÜMÜ

Bölüm 4.3'te SEY ile modelleme işleminde kullanılan şekil fonksiyonlarının tanımlanmasında, herbir elemanın düğüm numarası global düğüm numarası olarak ifade edilmiştir (Şekil 1). Karmaşık geometrili problemlerde şekil fonksiyonlarının global koordinat sisteminde verilmesi durumunda, hesaplama işleminin büyük sayılar üzerinde yapılması ve bununda hassasiyeti düşürmesinden dolayı, global koordinatlar yerine eleman üzerinde lokal koordinat sistemi tanımlanabilir. Şekil fonksiyonlarının lokal koordinat sistemi ifadeleri daha kısa ve hesaplamalar daha kolaydır.

Şekil 1'de bir üçgen eleman için verilen global koordinat sistemi ve bu elemana karşı gelen lokal koordinat sistemi gösterilmiştir. Burada, ( , )ξ η koordinat değişkenleri (0,1) aralığında değişir.

Şekil 1 Üçgen elemanda Global-Lokal koordinatlar: a. Üç düğüm noktalı bir üçgen elemanda global koordinatlar, b. Üç düğüm noktalı bir üçgen elemanda ( , )ξ η lokal koordinatlar ve değişim aralığı

Bir üçgen elemanın lokal koordinatlardaki şekil fonksiyonları;

1 2 3 1

φ ξ= φ η= φ = − − ξ η (3.1)

şeklinde ifade edilir. Şekil fonksiyonlarının toplamı bire eşittir.

123

1 2 1 1

1 .

i 0

i düğüm noktasında diğer düğüm noktalarında

φ φ φ

+ + =

= ⎨⎧

⎩

(3.2)

Lokal koordinatlar global koordinatlarla ilişkilendirilebilir.

1 1 2 2 3 3

x x x x

y y y y

φ φ φ

= + +

= + + (3.3)

13 23 3

x x x x

y y y y

ξ η

= + +

= + + (3.4)

( , )x y global ve ( , )ξ η lokal koordinat değerleri arasındaki dönüşüm dizey yapısında,

u x y u u

x x

u J u

u x y

y y

ξ ξ ξ

η η η

∂ ∂ ∂

⎧ ⎫ ⎡ ⎤ ∂⎧ ⎫ ⎧∂ ⎫

⎪∂ ⎪ ⎢∂ ∂ ⎥ ⎪ ⎪∂ ⎪∂ ⎪

⎪ ⎪=⎢ ⎥⎪ ⎪= ⎪ ⎪

⎨∂ ⎬ ⎢∂ ∂ ⎥⎨∂ ⎬ ⎨∂ ⎬

⎪ ⎪ ⎪ ⎪ ⎪ ⎪

⎢ ⎥ ∂ ∂

⎪∂ ⎪ ∂ ∂ ⎪⎩ ⎪⎭ ⎪⎩ ⎪⎭

⎩ ⎭ ⎣ ⎦

(3.5)

şeklindedir. Burada J Jacobian dizeyidir.

u u x

u u

ξ η

−

⎡∂ ⎤

⎢ ⎥

⎢∂ ⎥= ⎢∂ ⎥

⎢ ⎥

∂ ⎢∂ ⎥

⎢ ⎥

⎢∂ ⎥ ∂

⎣ ⎦ ⎣ ⎦

J (3.6)

(3.6) ifadesinde, A üçgen elemanın alanı olmak üzere, detJ = (x y_{13 23}−x y_{23 13}) 2= A ile ifade edilir.

124 EK 4 HibSis BİLGİSAYAR YAZILIMI

Bu bölümde ayrıntıları 6. bölümde verilen 2-boyutta elastik dalga yayılım benzeşimi gerçekleştiren "HibSis" isimli yazılımın GUI (Graphic User Interface) örnekleri verilmiştir.

Tez kapsamında geliştirilen HibSis, MATLAB sürüm 7.10.0.499 (R2010a) programlama dili kullanılarak kodlanmıştır (Şekil 1). HibSis yazılımı 2-Boyutta elastik dalga yayılımı modellemesi için geliştirilmiş olsa da yazılım ile sonlu fark, sonlu eleman ve melez model ağları oluşturulabilir. Kullanıcı arayüz aracılığıyla HibSis'in sunduğu tüm imkanlar kullanılabilir.

MATLAB programlama ortamının sunmuş olduğu hazır görselleştirme ve algoritma araçları nedeniyle HibSis yazılımının geliştirilmesinde tercih edilmiştir. Şekil 1' de sol panelde görüldüğü üzere HibSis üzerinde altı farklı işlem yapılabilir. Bunlar sırasıyla;

sonlu fark ağ oluşturma, sonlu eleman ağ oluşturma, ilk iki seçeneğe bağlı olarak ortak paylaşım bölgesi oluşturma, elastik dalga yayılımının sonlu fark benzeşimi, elastik dalga yayılımının sonlu eleman benzeşimi ve melez benzeşimidir.

125 Şekil 1 HibSis ana pencere

Şekil 2' de SFY’ne göre tanımlanan basit bir dörtgen ağ gösterilmiştir. Bu tür bir sonlu fark ağı Şekil 2'de sol panelde gösterilen giriş bilgileri kullanıcı tarafından tanımlanmak suretiyle kolayca yaratılabilir. Bu bilgiler sırasıyla; yatay (x) ve düşey (z) doğrultulardaki düğüm sayısı, yatay ve düşey doğrultulardaki uzamsal örnekleme aralıkları (m), başlangıç (referans) noktası yatay ve düşey koordinat değerleri ve sonlu fark ağ için bir model ad bilgilerinden oluşur.

126 Şekil 2 HibSis'de SFY'ne göre model ağı oluşturma

SFY ile oluşturulan ağ bilgileri gerekli ise yeniden düzenlenebilir. Şekil 3'de gösterildiği gibi, sağ panelde aktif sonlu fark ağına ilişkin model bilgileri yer alır.

Kullanıcı model bilgilerinde (örneğin, model boyutlarının, çok tabakalı bir ortamda her bir tabakaya ait sismik parametrelerin güncellenmesi gibi) güncelleme işlemi yaparak istenen geometri (düzenli geometri) tanımlayabilir (Şekil 3). Sonlu fark model ağ altında modele ilişkin gerekli renkli açıklama kullanıcı tarafından kolaylıkla takip edilebilir düzendedir. Kullanıcı bu noktada yeni bir modeli de kolaylıkla oluşturabilir.

Fakat bunu yapmak için bir önceki modelin temizlenmesi gerekir. Bu işlem "clear model" düğmesi ile kolaylıkla yapılabilir. Ayrıca “Select Source Location” düğmesine tıklanarak, kaynak yeri interaktif olarak istenilen şekilde belirlenebilir.

127

Şekil 3 HibSis'le hazırlanan sonlu fark model ağın, model koordinat ve eleman bilgilerinin düzenlenmesi.

Melez modelleme benzeşim işleminde altlık adı verilen sonlu fark ağı üzerinde istenilen takaka geometrisini oluşturabilmek için sonlu eleman yönteminden yararlanılır. Şekil 4' te bu durum gösterilmiştir. Sonlu fark ağ üzerine sonlu eleman yardımı ile fiziksel ve geometrik özellikleri farklı bir tabaka işlenmiştir. Bu şekilde model ağ bütününde oldukça karmaşık yapılar oluşturulabilir. Şekil 4'te görüleceği üzere kullanıcının kolaylıkla takip ve kontrol edebilmesi için sonlu eleman bilgileri eleman üzerinde gösterilmektedir. Burada global düğüm numaraları gösterilmiştir. Bu noktada, HibSis yazılımı sonlu fark ağ ile sonlu eleman ağ üzerindeki eleman ve düğümleri otomatik olarak tutar. Sonlu model ağ içerisinde karmaşık geometrili bir yapının sonlu eleman yöntemi ile oluşturulması fare yardımıyla oldukça kolaydır. Bunun için kullanıcı Şekil 4'te sarı renkle tanımlanan ağ bölgesi içerisinde sonlu eleman bölgesi sınırını fare

128

yardımıyla tanımlar. HibSis yazılımı sonlu model ağına ilişkin tüm tanımlamalarda SI birim sistemini kullanır.

Şekil 4 HibSis'te SEY'ne göre model ağı oluşturma

Sonlu farklar ağına benzer şekilde kullanıcı sonlu eleman ağı içinde gerekli gördüğünde güncelleştirme işlemi yapabilir. Şekil 5'te bu durum betimlenmiştir. Şekil 5'te görüldüğü gibi, model penceresi üzerinde sonlu eleman ağına ilişkin kullanılan bilgiler verilmektedir. Kullanıcı "start" düğmesine kliklediğinde modelin fiziksel ve geometrik özellikleri üzerinde gerekli güncelleme işlemini kolaylıkla yapılır. Kullanıcı bu noktada sonlu eleman tipini (dörtgen veya üçgen), sonlu ağ türünü (düzenli veya düzensiz), ağın yinelenip yinelenmiyeceğini, ağ giriş dosyalarının (koordinat ve eleman dosyaları) yolu ve adı, kaynak türünü (Ricker veya Gauss dalgacığı), ağ üzerinde kaynak lokasyonunu, yineleme sayısını, lokal-global koordinat dönüşüm işleminin yapılıp yapılmayacağını ve

129

çözücü (system solver) türünü değiştirebilir. Şekil 5'te sonlu eleman ağ modeli üzerinde global düğüm noktaları ile birlikte sismik hız (Vp) renk ölçeğinde verilmektedir.

Şekil 5 HibSis'te sonlu eleman model ağında model bilgilerinin düzenlenmesi

2-Boyutta elastik dalga yayılımının melez yöntemle benzeşimi işleminin son adımı, sonlu fark ve sonlu eleman ağ için bir ortak bölge oluşturulmasıdır. Ortak bölge, sonlu fark ve sonlu eleman ağ tarafından kullanılan ortak düğümlerin değerlerinin her yineleme sonunda değiş tokuşunun yapıldığı alanı tanımlar. Kullanıcı ortak bölge için iki farklı yol izleyebilir. Birincisi ortak bölgenin dikdörtgen geometride oluşturulması, bir diğeri ise poligon geometride oluşturulmasıdır. Ortak bölgenin her iki geometride oluşturulması tümüyle fare yardımıyla gerçekleştirilir. Şekil 6'da örnek bir ortak bölge sonlu elemanlara ait global düğüm noktalarıyla gösterilmiştir. Geliştirilen HibSis

130

yazılımı ile 2-Boyutta elastik dalga yayılımı benzeşimi melez yöntem kullanılarak yapılabileceği gibi yalnızca SFY veya yalnızaca SEY kullanılarakta yapılabilir.

Şekil 6 HibSis'te ortak bölgenin oluşturulması

Sarı bölge SFY uygulanan alan, pembe bölge SEY uygulanan alan, mavi bölge ortak bölge (Transition Zone)

HibSis Programı Model Dosya Örnekleri

HibSis programı sonlu farklar ağının oluşturulması için üç adet dosya kullanır. Bu dosyalardan ilki bilgi dosyasıdır ve modele ilişkin genel tanımlamaları içerir (Çizelge 1). İkinci dosya model eleman bilgilerinin tutulduğu dosyadır (Çizelge 2). Üçüncü dosya ise model düğüm bilgilerinin tutulduğu dosyadır (Çizelge 3). Bu dosyaların içerikleri ve ilgili dosyaları okuyan MATLAB fonksiyonları aşağıdaki verilen yapıdadır.

131 Çizelge 1 Bilgi dosyası içeriği

basamakh.dat

model_header 2401 2500 1 150 0 300 10 10 50.0 1000 2000 1 1

basamake.dat basamakn.dat

satır-1 Model isimi

satır-2 Modelde kullanılan düğüm sayısı, eleman sayısı ve yoğunluk

satır-3 eleman numarası, başlangıç zamanı (ms), bitiş zamanı (ms), kaynağın x koordinatı, kaynağın z koordinatı, kaynak frekansı (Hz)

satır-4 minumum ve maxsimum hız, kaynak için katsayı, model için katsayı

satır-5 eleman bilgi dosya ismi satır-6 düğüm bilgi dosya ismi

Çizelge 2 Eleman bilgi dosyası içeriği basamake.dat

model_eleman

en nd1 nd2 nd3 nd4 rho Vp Vs 1 1 51 52 2 2100.00 2000.00 1000.00 2 2 52 53 3 2100.00 2000.00 1000.00 3 3 53 54 4 2100.00 2000.00 1000.00 4 4 54 55 5 2100.00 2000.00 1000.00 5 5 55 56 6 2100.00 2000.00 1000.00 .

. .

2398 2446 2496 2497 2447 2700.00 4000.00 2500.00 2399 2447 2497 2498 2448 2700.00 4000.00 2500.00 2400 2448 2498 2499 2449 2700.00 4000.00 2500.00 2401 2449 2499 2500 2450 2700.00 4000.00 2500.00 satır-1 Model isimi

satır-2 eleman numarası, düğüm numaraları, yoğunluk ve hızlar sütun-1 Eleman numarası

sütun-2 Eleman birinci düğüm numarası sütun-3 Eleman ikinci düğüm numarası sütun-4 Eleman üçüncü düğüm numarası sütun-5 Eleman dördüncü düğüm numarası sütun-6 Eleman yoğunluk değeri (gr/cm³)

sütun-7 Eleman sismik boyuna dalga (P) hızı (Vp) (m/s) sütun-8 Eleman sismik makaslama dalga (S) hızı (Vs) (m/s)

132 Çizelge 3 Düğüm bilgi dosyası içeriği

basamakn.dat model_node

n x z vx0 vz0 cnd vh 1 0.0 0.0 0 0 0 0 2 20.0 0.0 0 0 0 0 3 40.0 0.0 0 0 0 0 4 60.0 0.0 0 0 0 0 5 80.0 0.0 0 0 0 0 6 100.0 0.0 0 0 0 0 7 120.0 0.0 0 0 0 0 8 140.0 0.0 0 0 0 0 .

. .

2498 940.0 980.0 0 0 0 0 2499 960.0 980.0 0 0 0 0 2500 980.0 980.0 0 0 0 0 satır-1 Model isimi

satır-2 eleman numarası, x ve z koordinatları sütun-1 Düğüm numarası

sütun-2 Düğüm X koordinatı sütun-3 Düğüm Z koordinatı

sütun-4 Düğüm noktasının başlangıç yatay değeri sütun-5 Düğüm noktasının başlangıç düşey değeri

sütun-6 Düğüm noktasına ait sınır koşulu (0= serbest düğüm, 1= yatay doğrultuda kuvvet uygulanır, düşey yönde düğüm serbest, 2= düşey doğrultuda kuvvet uygulanır, yatay yönde düğüm serbest, 3= yatay ve düşey yönde kuvvet uygulanır)

sütun-7 Düğüm noktasının geçmiş değerleri (0=tutulmasın,1=tutulsun)

133

% SFY için bilgi dosyasını, model eleman ve düğüm dosyalarını okuyan fonksiyon function [mtitle,nel,ndpt,n,st,et,sx,sz,f,vmin,vmax,scoef,mcoef,grav,filename]=altlik [filein,mfilepath] = uigetfile({'*.DAT','DAT Files [*.DAT]';

'*.TXT','TXT Files [*.TXT]';

'*.*' ,'All Files [*.*]'},'Input Model File') ; [pathstr,fname,ext,ver] = fileparts(filein);

fid = fopen(filein,'rt') ;

mtitle = fgetl(fid) ; % model baslıgını oku tmp = fscanf(fid,'%u%u%f',[1,3]) ; % pass gec

nel = tmp(1,1); ndpt = tmp(1,2); % eleman ve toplam düğüm sayısı grav = tmp(1,3);

tmp=fscanf(fid,'%f%f%f%f%f%f',[1,6]);

n = tmp(1,1) ; st = tmp(1,2) ; % eleman numarası,baslangıç zamanı(ms) et = tmp(1,3); sx=tmp(1,4); % bitiş zamanı(ms), kaynagın x koordinatı sz=tmp(1,5); % kaynagın z koordinatı, kaynak frekansı(Hz) f = tmp(1,6);

tmp = fscanf(fid,'%f',[1,4]) ;

vmin = tmp(1,1); vmax = tmp(1,2); % minumum ve maxsimum hız

scoef = tmp(1,3); mcoef = tmp(1,4); % kaynak için katsayı, model için katsayı ename = fscanf(fid,'%s',[1,1]) ; % eleman bilgi dosya ismi,

nname = fscanf(fid,'%s',[1,1]) ; % node bilgi dosya ismi fclose(fid) ;

clear fid tmp ans mstruct ;

import_eleman(ename);

import_node(nname);

filename=[fname 'fem' ext];

end

function import_eleman(ename)

% SFY için eleman dosyasını okur global el

eofstate=0;

fid=fopen(ename,'rt');

hdr=fgetl(fid);

tline=fgetl(fid);

while ~eofstate

tmp = fscanf(fid,'%u%u%u%u%u%f%f%f',[1,8]);

n = tmp(1) ;

el(n).node = [tmp(2),tmp(3),tmp(4),tmp(5)] ; el(n).rho = tmp(6);

el(n).vhv = [tmp(1,7), tmp(1,8)];

eofstate = feof(fid) ; end;

fclose(fid);

end

134 function import_node(nname)

% SFY için düğüm dosyasını okur global node

eofstate = 0 ;

fid = fopen(nname,'rt') ; mhdr = fgetl(fid);

tline = fgetl(fid);

while ~eofstate

tmp = fscanf(fid,'%u%f%f%u%u%u%u',[1,7]) ; n = tmp(1) ;

node(n).coord = [tmp(2), tmp(3)] ; node(n).vhv0 = [tmp(4), tmp(5)] ; node(n).nc = tmp(6);

node(n).vhist= tmp(7);

eofstate = feof(fid) ; end;

fclose(fid) ; end

HibSis programı, sonlu eleman ağının oluşturulması için iki adet dosya kullanır. Bu dosyalardan ilki, model eleman bilgilerinin tutulduğu dosyadır (Çizelge 4). İkincisi ise model düğüm bilgilerinin tutulduğu dosyadır (Çizelge 5). Bu dosyaların içerikleri ve ilgili dosyaları okuyan MATLAB fonksiyonları aşağıdaki verilen yapıdadır.

135 Çizelge 4 Eleman bilgi dosyası içeriği

eleman_new.dat

en nd1 nd2 nd3 rho Vp Vs 1 21 3 246 2700 4000 2500 2 17 2 166 2100 2000 1000 3 16 10 293 2700 4000 2500 4 11 1 29 2100 2000 1000 5 66 11 102 2700 4000 2500 6 1 5 29 2100 2000 1000 7 11 12 102 2700 4000 2500 8 102 12 290 2700 4000 2500 . . .

satır-1 eleman numarası, düğüm numaraları, yoğunluk ve hızlar sütun-1 Eleman numarası

sütun-2 Eleman birinci düğüm numarası sütun-3 Eleman ikinci düğüm numarası sütun-4 Eleman üçüncü düğüm numarası sütun-5 Eleman yoğunluk değeri (gr/cm³)

sütun-6 Eleman sismik boyuna dalga (P) hızı (Vp) (m/s) sütun-7 Eleman sismik makaslama dalga (S) hızı (Vs) (m/s)

Çizelge 5 Düğüm bilgi dosyası içeriği coord_new.dat

nn X Z Vp Vs 1 0.00000 260.00000 2000 1000 2 520.00000 260.00000 2000 1000 3 520.00000 540.00000 4000 2500 4 980.00000 540.00000 4000 2500 5 0.00000 220.00000 4000 2500 6 980.00000 220.00000 4000 2500 7 980.00000 980.00000 4000 2500 8 0.00000 980.00000 4000 2500 9 490.00000 220.00000 2000 1000 . . .

satır-1 Düğüm numarası, X ve Z koordinatları ve hızlar sütun-1 Düğüm numarası

sütun-2 Düğüm X koordinatı

sütun-4 Düğüm sismik boyuna dalga (P) hızı (Vp) (m/s) sütun-5 Düğüm sismik makaslama dalga (S) hızı (Vs) (m/s)

136 function fem_model

% SEY için model ağı oluşturan fonksiyon

[mtitle,nel,ndpt,n,st,et,sx,sz,f,vmin,vmax,scoef,mcoef,grav,filename]=altlik;

XY=[2,2;4,6;0,0;0,520;980,520;980,980;0,980;980,0;980,260;0,260;0,540;0,540;0,980;

0,980;];

d = decsg(XY);

[p,e,t] = initmesh(d);

[p,e,t] = refinemesh(d,p,e,t);

v = zeros(length(t));

i1 = pdesdt(t,1);

vp(i1) = 2000;

vs(i1) = 1000;

ro(i1) = 2100;

i2 = pdesdt(t,2);

vp(i2) = 4000;

vs(i2) = 2500;

ro(i2) = 2700;

createMesh_coord_new(coord,eleman,vp,vs);

hold on

pdeplot(p,e,t,'xydata',eleman_new(:,6),'colorbar','on','xystyle','flat','mesh','on');

axis ij

for i=1:length(p(1,:)) hold on

text(p(1,i),p(2,i),num2str(i),'FontSize',10,'fontweight','bold') end

end

HibSis programı melez (Hibrid) modelleme ağının oluşturulması için Çizelge 1-3’te verilen SFY model bilgi, eleman ve düğüm dosyaları ve Çizelge 4-5’te verilen SEY eleman ve düğüm dosyalarını kullanılır. Her iki yöntemin ortak bölgedeki (transition zone) ortak düğümlerinin tutulduğu dosya örneği Çizelge 6’daki gibidir. Bu dosyaların içerikleri ve ilgili dosyaları okuyan MATLAB fonksiyonu aşağıdaki verilen yapıdadır.

137

Çizelge 6 Melez modelleme için SFY ve SEY ortak düğüm ve eleman bilgi dosyası içeriği

FDFE_ortak.dat

FD_nd ortFE_nd xFE zFE FE_el_num FE_nd1 FE_nd2 FE_nd3 551 5 0.00 220.00 6 1 5 29 552 29 28.99 220.00 4 11 1 29 552 29 28.99 220.00 6 1 5 29 552 29 28.99 220.00 29 11 29 30 554 30 66.67 220.00 9 12 11 30 554 30 66.67 220.00 11 13 12 30 554 30 66.67 220.00 29 11 29 30 554 30 66.67 220.00 31 13 30 31 557 31 126.11 220.00 26 14 13 31 557 31 126.11 220.00 31 13 30 31 557 31 126.11 220.00 33 14 31 32 . . .

satır-1 Geçiş bölgesindeki SFY düğüm numarası, SEY düğüm numarası, ilgili SEY düğüm numarasının X ve Z koordinatları, ilgili SEY düğümünün eleman numarası ve bu elemanın düğüm numaraları sütun-1 SFY düğüm numarası

sütun-2 SEY düğüm numarası

sütun-3 ilgili SEY düğüm numarasının X koordinatı sütun-4 ilgili SEY düğüm numarasının Z koordinatı sütun-5 ilgili SEY düğümünün eleman numarası sütun-6 elemanın birinci düğüm numarası sütun-7 elemanın ikinci düğüm numarası sütun-8 elemanın üçüncü düğüm numarası

138 function createMesh_FDFE

% SFY ve SEY model dosyalarını kullanarak Melez model ağı oluşturur

% Sonlu Farklar için model düğüm ve elemanlarının belirlendiği kısım createMesh_altlik(nx,nz,dx,dz,x0,z0,fname)

[mtitle,nel,ndpt,n,st,et,sx,sz,f,vmin,vmax,scoef,mcoef,grav,filename]=altlik;

drawmodel(nel);

% Sonlu Elemanlar için model düğüm ve elemanlarının belirlendiği kısım hold on

NumberOfSub=1;

% basamak fay

XY = [2;6;0;520;520;980;980;0;260;260;540;540;980;980;];

plot(XY(3),XY(9),'*r',XY(4),XY(10),'*r',XY(5),XY(11),'*r',XY(6),XY(12),'*r',XY(7), XY(13),'*r',XY(8),XY(14),'*r','LineWidth',3);

% Sonlu Farklar ve Sonlu Elemanlar için ortak düğümlerin belirlendiği kısım

% Geçiş-bölgesi (Average_Zone) AveZones=struct('X',[0 0 0 0 ], ...

'Y',[0 0 0 0 ]);

AverageZone;

yx1 = AveZones.X(1,1); yx2 = AveZones.X(1,2); yx3 = AveZones.X(1,3); yx4 = AveZones.X(1,4);

yz1 = AveZones.Y(1,1); yz2 = AveZones.Y(1,2); yz3 = AveZones.Y(1,3); yz4 = AveZones.Y(1,4);

XY = [3; 4; yx1; yx2; yx3; yx4; yz1; yz2; yz3; yz4];

d = decsg(XY);

[p,e,t] = initmesh(d);

query = questdlg('Alt bölge daha sık üçgenlensin mi ?', ...

'FEM domain ','YES', 'NO', 'NO');

switch (upper(query)) case 'YES'

ic_sub = find(t(4,:)==2);

tt = t(:,ic_sub);

[p,e,t]=refinemesh(d,p,e,tt);

end

v = zeros(length(t));

i1 = pdesdt(t,1);

vp(i1) = 2000;

vs(i1) = 1000;

ro(i1) = 2100;

i2 = pdesdt(t,2);

vp(i2) = 4000;

vs(i2) = 2500;

ro(i2) = 2700;

q=pdetriq(p,t);

pdeplot(p,e,t,'xydata',vp,'colorbar','on','xystyle','flat','mesh','on');

for i=1:length(p(1,:)) hold on

139

text(p(1,i),p(2,i),num2str(i),'FontSize',10,'fontweight','bold') end

% Average zone daki Sonlu Farklar düğümlerinin Sonlu Elemanlar daki

% hangi elemana denk geldiğinin bulunması k=[];

r=zeros(1,length(p));

rm=1;

fid2 = fopen('FDFE_ortak.dat','wt');

fprintf(fid2,'%s\n','FD_nd ortFE_nd xFE zFE FE_el_num FE_nd1 FE_nd2 FE_nd3');

for ic = AveNodes.Node(1): AveNodes.Node(3) for j = 1:length(p)

if (( AveZones2.X(ic)+ dx/2 ) > p(1,j) ) && (( AveZones2.X(ic)- dx/2 ) < p(1,j) )

&& ...

(( AveZones2.Y(ic)+ dz/2 ) > p(2,j) ) && (( AveZones2.Y(ic)- dz/2 ) < p(2,j) ) k=[ k j];

end end

r= sqrt((AveZones2.X(ic)-p(1,k)).^2 + (AveZones2.Y(ic)-p(2,k)).^2);

[v,ind]=min(r);

j=k(ind);

if isempty(j), continue; end for i=1:length(t)

for L=1:3

if (find(t(L,i)==j))

fprintf(fid2,'%u\t %u\t %4.2f\t %4.2f\t %u\t %u\t %u\t %u\t\n',ic, j, p(1,j), p(2,j),i,t(1,i),t(2,i),t(3,i));

fdd(rm,1) = ic;

fdd(rm,2) = i;

rm = rm+1;

end end end k=[];

end

fclose(fid2);

[FD_nd,ortFE_nd,xFE,zFE,FE_el_num,FE_nd1,FE_nd2,FE_nd3]=textread('FDFE_orta k.dat','%d%d%f%f%d%d%d%d','headerlines',1);

% Sonlu Farklar ve Sonlu Elemanlar bilgi dosyası fid3 = fopen('FD-FE_bilgi_dosyası.dat','wt');

fprintf(fid3,'%s\n','FD için hazırlanna altlık ve üzerine eklenene FE ile ilgili bilgiler ');

fprintf(fid3,'%u\t\t\t %s\n',nx,'nx number of nodes in x-direction ');

fprintf(fid3,'%u\t\t\t %s\n',nz,'nz number of nodes in z-direction ');

fprintf(fid3,'%u\t\t\t %s\n',dx,'dx sampling rate in x-direction ');

fprintf(fid3,'%u\t\t\t %s\n',dz,'dz sampling rate in z-direction ');

fprintf(fid3,'%u\t\t\t %s\n',x0,'x0 centre location in x-direction ');

140

fprintf(fid3,'%u\t\t\t %s\n',z0,'z0 entre location in z-direction ');

fprintf(fid3,'%s\t\t%s\n',fname, 'model file name ');

fname1 = [fname,'h.dat'];

fname2 = [fname,'n.dat'];

fname3 = [fname,'e.dat'];

fname4 = [fname,'_FD.dat'];

fprintf(fid3,'%s\t\t%s\n',fname1, 'Sonlu Farklar header dosyası ' );

fprintf(fid3,'%s\t\t%s\n',fname2, 'Sonlu Farklar node dosyası ' );

fprintf(fid3,'%s\t\t%s\n',fname3, 'Sonlu Farklar eleman dosyası ' );

fprintf(fid3,'%s\t\t%s\n',fname4, 'Sonlu Farklar hesaplamada kullanılacak Hız içeren dosyası ' );

fprintf(fid3,'%u\t\t\t %s\n',nel, 'nel FD number of eleman ');

fprintf(fid3,'%u\t\t\t %s\n',ndpt, 'ndpt FD number of node ');

fprintf(fid3,'%u\t\t\t %s\n',grav, 'grav FD density ');

fprintf(fid3,'%u\t\t\t %s\n',n, 'n FD iteration ');

fprintf(fid3,'%u\t\t\t %s\n',st, 'st FD start time ');

fprintf(fid3,'%u\t\t\t %s\n',et, 'et FD end time ');

fprintf(fid3,'%u\t\t\t %s\n',sx, 'sx FD source x-cordinate');

fprintf(fid3,'%u\t\t\t %s\n',sz, 'sz FD source z-cordinate');

fprintf(fid3,'%u\t\t\t %s\n',f, 'f FD source frequency ');

fprintf(fid3,'%u\t\t\t %s\n',vmin, 'vmin FD minimum velocity ');

fprintf(fid3,'%u\t\t\t %s\n',vmax, 'vmax FD maximum velocity ');

fprintf(fid3,'%s\t\t\t %s\n','eleman.dat', 'Sonlu Elemanlar eleman dosyası ');

fprintf(fid3,'%s\t\t\t %s\n','coord2.dat', 'Sonlu Elemanlar node dosyası ');

fprintf(fid3,'%u\t\t\t %s\n',FE_nel, 'Sonlu Elemanlar number of nodes');

fprintf(fid3,'%u\t\t\t %s\n',FE_ndpt, 'Sonlu Elemanlar number of eleman');

fprintf(fid3,'%s\t %s\n','FDFE_boundary.dat', 'Average Zone daki sınır node lar ');

fprintf(fid3,'%s\t %s\n','FDFE_ortak.dat', 'FD,FE ve average zone için ortak node ve elemanlar');

fclose(fid3);

end

141 Çizelge 7 Melez modelleme bilgi dosyası

FD için hazırlanna altlık ve üzerine eklenene FE ile ilgili bilgiler 50 nx number of nodes in x-direction 50 nz number of nodes in z-direction 20 dx sampling rate in x-direction 20 dz sampling rate in z-direction 0 x0 centre location in x-direction 0 z0 entre location in z-direction basamak model file name

basamakh.dat Sonlu Farklar header dosyası basamakn.dat Sonlu Farklar node dosyası basamake.dat Sonlu Farklar eleman dosyası

basamak_FD.dat Sonlu Farklar hesaplamada kullanılacak Hız içeren dosyası 2401 nel FD number of eleman

2500 ndpt FD number of node 1 grav FD density 300 n FD iteration 0 st FD start time 300 et FD end time 5 sx FD source x-cordinate 5 sz FD source z-cordinate 50 f FD source frequency 1000 vmin FD minimum velocity 4000 vmax FD maximum velocity eleman.dat Sonlu Elemanlar eleman dosyası coord.dat Sonlu Elemanlar node dosyası 390 Sonlu Elemanlar number of nodes 713 Sonlu Elemanlar number of eleman FDFE_boundary.dat Average Zone daki sınır node lar

FDFE_ortak.dat FD,FE ve average zone için ortak node ve elemanlar

142

% Ricker Kaynak Fonksiyonu ts = 0.02; % fs = 50 Hz dt = 0.001;

t = (1 : 1024)*dt;

ls = 3*ts/dt;

tg = pi*(t(1 : ls)-1.5*ts) / (1.5*ts);

src = (1-4*tg.^2).*exp(-2*tg.^2);

src((ls+1) : 1024) = 0*((ls+1) : 1024);

t = t*1000;

plot (t(1 : ls), src(1 : ls));

title('Ricker '); xlabel('Zaman [ms]'); ylabel('Genlik');

Şekil 7 Ricker kaynak fonksiyonu

Belgede ANKARA ÜNİVERSİTESİ FEN BİLİMLERİ ENSTİTÜSÜ DOKTORA TEZİ MELEZ ( HİBRİD ) YÖNTEMLE ELASTİK DALGA DENKLEMİNİN İKİ-BOYUTTA MODELLENMESİ (sayfa 118-154)