エカつきのブログ = Eka tsuki no blog

Program Matlab Persamaan Transport 2 Dimensi jika diketahui koefisien afeksi (laju) dan Koefisien difusi pada arah x dan y

clc;
clear;
disp('----------------------------------------------------')
disp('Program Persamaan Transport 2 Dimensi non-dimesional')
disp(' dengan koefisien afeksi(laju) dan koefisien difusi')
disp('----------------------------------------------------')
disp('syarat awal  : f(x)=sin(x)*cos(y)')
disp('syarat batas : u(x,0,t)= 0')
disp('               u(x,1,t)= 0')
disp('               u(0,y,t)= 0')
disp('               u(1,y,t)= 1')
disp(' ')
%% Definisi variabel Non-Dimensional

X=1; %panjang arah x 
Y=1; %lebar arah y
T=1; %waktu akhir

%% Input Nilai Konstanta dan Banyaknya Segmen

disp('Input nilai koefisien afeksi dan koefisien difusi')
flag=0;
while flag==0
    alpha_x=input('Masukkan koefisien afeksi(laju) alpha(x) = ');
    if alpha_x<=0
        disp('Koefisien alpha(x) harus positif')
    else
        break
    end
end

while flag==0
    alpha_y=input('Masukkan koefisien afeksi(laju) alpha(y) = ');
    if alpha_y<=0
        disp('Koefisien alpha(y) harus positif')
    else
        break
    end
end

while flag==0
    beta_x=input('Masukkan koefisien difusi beta(x) = ');
    if beta_x<=0
        disp('Koefisien beta(x) harus positif')
    else
        break
    end
end

while flag==0
    beta_y=input('Masukkan koefisien difusi beta(y) = ');
    if beta_y<=0
        disp('Koefisien beta(y) harus positif')
    else
        break
    end
end

disp(' ')
disp('Input banyaknya segmen')
while flag==0
    k=input('Masukkan banyaknya segmen pada x dan y = ');
    if k<=0 || ceil(k)-k>0
        disp('Banyaknya segmen harus bulat positif')
    else
        break
    end
end

while flag==0
    maxn=input('Masukkan banyaknya segmen waktu = ');
    if maxn<=0 || ceil(maxn)-maxn>0
        disp('Banyaknya segmen waktu harus bulat positif')
    else
        break
    end
end

dx=X/k;
dy=Y/k;
dt=T/maxn;

%% Syarat Kestabilan

P=beta_x*dt/dx^2;
Q=beta_y*dt/dy^2;
R=alpha_x*dt/dx;
S=alpha_y*dt/dy;
disp(' ')
if P+Q+R+S<=1/2
    disp('Kondisi Kestabilan : Stabil')
else
    disp('Kondisi Kestabilan : Tidak stabil')
end

%% Grid Komputasi

for j=1:k+1
    y(j)=(j-1)*dy;
    for i=1:k+1
        x(i)=(i-1)*dx;
        u(i,j,1)=sin(pi*x(i))*cos(pi*y(j)); %syarat awal
    end
end

for n=1:maxn+1;
    waktu(n)=(n-1)*dt;
    for j=1:k+1
        for i=1:k+1
            u(1,j,n)=0; %syarat batas
            u(i,1,n)=0;
            u(i,k+1,n)=0;
            u(k+1,j,n)=1;
        end
    end
end

%% Penyelesaian dengan metode FTCS

for n=1:maxn
    for j=2:k
        for i=2:k
            u(i,j,n+1)=(P+R)*u(i+1,j,n)+(Q+S)*u(i,j+1,n)+P*u(i-1,j,n)+Q*u(i,j-1,n)+(1-2*P-2*Q-R-S)*u(i,j,n);
        end
    end
end

%% output plot tiga dimensi keadaan panas pada akhir waktu

mesh(x,y,u(:,:,maxn)')
xlabel('Panjang')
ylabel('Lebar')
zlabel('Suhu') 
title('Grafik Penyebaran Panas Persamaan Transport 2 Dimensi')