Simulasi Numerik,NUMERICAL MODELLING : GELOMBANG AKUSTIK (Tanpa batas penyerap)
NUMERICAL MODELLING : GELOMBANG AKUSTIK
Di sini kita akan
menggunakan metode numerik centered
difference staggered grid.
Persamaan
di atas disebut dengan centered finite
difference time domain 2D orde satu dengan akurasi spasial orde dua. Lalu, persamaan
penjalaran gelombang akustik untuk medium 2D yang akan kita selesaikan dengan
metode numeris di atas:
P di
atas merupakan pressure, vx dan vz
merupakan kecepatan partikel ke arah x dan z. ρ merupakan densitas medium dan
ca merupakan kecepatan gelombang P. Nilai dx, dz, dan dt merupakan
sampling spasial arah x, z, dan sampling waktu. Kemudian
persamaan (2) diselesaikan menggunakan pers (1) dengan memisalkan medium
kita berukuran Z*X, menjadi :
P(i,j) = P(i,j)-dt*ρ(i,j)*ca(i,j)^2*((vx(i+1,j)-vx(i,j))/dx+(vz(i,j+1)-vz(i,j))/dy)
vx(i,j)= vz(i,j)-dt*(P(i,j)-P(i-1,j))/dx/ρ(i,j)
vz(i,j)= vx(i,j)-dt*(P(i,j)-P(i,j-1))/dy/ρ(i,j)
(i=1,2….,Z; j=1,2,…,X.) (3)
Dari
persamaan di atas kita ketahui bahwa untuk menyelesaikan P, vx dan vz, kita
membutuhkan parameter rho dan kecepatan P (ca) terlebih dahulu.
Membuat Parameter
Model untuk simulasi
Untuk membuat matriks parameter
model lapisan datar tentu sangat gampang, cukup dengan perulangan menggunakan
perintah looping seperti
for..end.
misal :
Nx = 201;
Ny = 101;
for i=1:Nx
for j=1:floor(Ny/2);
rho(i,j) = 1800;
c(i,j) = 2000;
end
for j=floor(Ny/2):Ny;
rho(i,j) = 2300;
c(i,j) = 2600;
K(i,j)=rho(i,j)*c^2;
end
end
Bagaimana dengan model seperti ini?
Tentu susah bukan jika menggunakan perulangan biasa? untuk
itu kita perlu akali dengan cara lain seperti berikut :
1.
Pertama, gambar sendiri model menggunakan image
editor software seperti Paint, photoshop, atau Paint.NET (recommended).
2.
Set ukurannya sesuai model kita, misal untuk
tugas kita disuruh membuat model 2000 x 1000 m dengan jarak spasial dx/dz=10 m.
Maka ukuran matriksnya (ukuran model dibagi jarak spasial) sebesar 200x100.
dan matikan anti-aliasing.
3.
Kemudian gambar model anda dengan warna
grayscale (warna dimana nilai R, G, dan B-nya sama dari 0 sampai 255).
misal lapisan 1 R: 233 B:233 G:233, lapisan
2 : R: 213 B:213 G:213.
Kemudian save dengan format BMP dengan
setting bit-depth = 8 bit.
4.
Kemudian buka di matlab dengan perintah :
data=imread('namafile.bmp');
Sekarang variabel ‘data’ sudah terisi tanda
yang menunjukkan masing-masing lapisan, misal untuk model di atas, lapisan satu
pada matriks ditandai dengan nilai 0, sedangkan lapisan 2 bernilai 1.
Untuk pembuatan model, tinggal dikondisikan
saja sesuai ketentuan :
misal lapisan 1 vp=2300 &
rho=2000, lapisan 2 vp=2500 & rho=2200
Nx = 200;
Nz = 100;
for i=1:Nz
for j=1:Nx
if data(i,j)==0
rho(i,j) = 2000;
c(i,j) = 2300;
else if data(i,j)==1
rho(i,j) = 2200;
c(i,j) = 2500;
end
end
end
end
lalu simpan data tersebut. Salah satu
contohnya dengan menggunakan perintah
>>
dlmwrite('c.txt',c,'delimiter','\t')
>> dlmwrite('rho.txt',rho,'delimiter','\t')
Program utama
simulasi gelombang seismik
Kita awali dengan membuat m-file baru, setelah itu seperti
biasa kita deklarasikan penghapusan memori dan penutupan figure-figure.
close all;
clear all;
clc
Setelah itu, kita deklarasikan parameter-parameter awal :
%input
Nx = 200; %panjang grid model
Ny = 100; %kedalaman grid model
nstep=300; %jumlah iterasi
f0=20; %frekuensi source dominan ricker
dx = 10; %jarak spasial horizontal
dy = 10; %jarak spasial vertikal
c=dlmread('c.txt'); %matriks vp model
rho=dlmread('rho.txt'); %matriks rho model
Lalu,
kita buat posisi koordinat untuk penampilan simulasi dan penghitungan parameter.
Parameter yang diload tadi kita transpose. Kenapa dilakukan transpose? format
pengindeksan matlab untuk array misal
a(i,j) dimana i merupakan menunjukkan posisi indeks vertikal kadang sering
tertukar oleh penulisan rumus yang terbiasa menulis i pada koordinat bagian
kiri sebagai posisi indeks horizontal.
%penghitungan parameter model dan posisi untuk koordinat
surf
c=c';
rho=rho';
K=rho.*c.^2;
is=1;
x1 = (0:Nx-2)*dx+dx/2;
y1 = (0:Ny-2)*dy+dy/2;
dt=1./sqrt(1./dx^2+1./dy^2)/max(max(c));
dt
Kemudian, kita pesan
memori awal dan menghitung koordinat untuk tampilan surf simulasi kita.
%pemesanan memori awal
p=zeros(Nx-1,Ny-1); %pressure
u=zeros(Nx,Ny-1); %vx
v=zeros(Nx-1,Ny); %vz
x=zeros(Nx-1,Ny-1);
y=zeros(Nx-1,Ny-1);
rec=zeros(nstep,Nx);
for j=1:Ny-1
x(:,j)=x1(:);
end
for i=1:Nx-1
y(i,:)=y1(:);
end
surf(x,y,p); view(2); axis equal; shading interp;
time=0.0;
Setelah
itu, perhitungan utama simulasi kita muat dalam program. Tiap perhitungan
matriks P, u, dan v, iterasi ditambahkan gangguan berupa source ricker wavelet dan hasil perhitungan pada elemen matriks
komponen v pada posisi permukaan model disimpan ke dalam variabel rec sebagai receiver. Untuk tampilan
animasi, tiap 10 iterasi kita kondisikan agar data hasil terakhir perhitungan
kita surf.
for time_steps = 1: nstep
for j=1:Ny-1
for i=1:Nx-1;
p(i,j)=p(i,j)-dt*K(i,j)*((u(i+1,j)-u(i,j))/dx+(v(i,j+1)-v(i,j))/dy);
end
end
for j=1:Ny-1
for i=2:Nx-1;
u(i,j)=u(i,j)-dt*(p(i,j)-p(i-1,j))/dx/rho(i,j);
end
end
for j=2:Ny-1
for i=1:Nx-1;
v(i,j)=v(i,j)-dt*(p(i,j)-p(i,j-1))/dy/rho(i,j);
end
end
%sumber
seismik
a=pi^2*f0^2;
f=1e8*2.0*a*(time-1.2/f0)*exp(-a*(time-1.2/f0)^2); %ricker
wavelet
v(is,Ny)=v(is,Ny)+dt*f/rho(is,1); %source diletakkan di pojok kiri
model
rec(time_steps,:)=u(:,Ny-2); %receiver diletakkan di permukaan
model
time=time+dt;
if
(mod(time_steps,10) == 0)
figure(gcf)
surf(x,y,p); view(2);
axis equal; shading interp;
fprintf(2,' time = %8.5f\n', time);
pause(0.1)
end
end
Setelah program tersebut dijalankan, rekaman seismik dapat
ditampilkan dengan perintah :
clim=[-500,500];
imagesc(x1,0:dt:nstep*dt,rec,clim);
colormap(flipud(gray));
ylabel('waktu (s)'); xlabel('jarak (m)');
Sekian :)
tips : karena tidak diberi kondisi batas penyerap, sebaiknya
iterasi (pada program atas disebut nstep) dikurangi jangan sampai gelombang
langsung pada permukaan memantul kembali dari bagian kiri model.
0 comments: