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:

(2)

P di atas merupakan pressure, vx dan vz merupakan kecepatan partikel ke arah x dan z. ρ merupakan densitas medium dan c_a 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 :

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

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.