50 Lampiran 1: Iterasi Newton
1
format long % agar memberikan nilai sebanyak 15 bilangan decimal
2 3
% nilai perkiraan awal
4
x= xguess;
5
y= yguess;
6 7
% input jumlah iterasi yang diinginkan
8
Niter=input('jumlah iterasi= ');
9
iterasiold(1) = x;
10
iterasiold(2) = y;
11
for n=1:Niter
12
f1=xˆ2-yˆ2-y;
13
f2=yˆ2+xˆ2-x;
14
j=[2*x-2*y-1; 2*x-1 2*y]; iterasi = [x; y]-inv(j)*[f1; f2]
15
galat = norm(iterasi-iterasiold')
16 17
x = iterasi(1);
18
y = iterasi(2);
19
iterasiold = iterasi';
20
end
Universitas Sumatera Utara
51 Lampiran 2: Aliran Satu Fase
1
maxk = 10000;
2
tfinal = 200;
3
dt = tfinal/maxk;
4
n = 30;
5
dx = 1/n;
6
b=dt/2*dx;
% jumlah waktu
% jumlah diskritisasi
7 8
% Initial condition
9
x(1) = 0.0;
10
u(1,1) = 1.0;
11 12
for i = 2:n+1
13
x(i) =x(i-1) + dx;
14
u(i,1) =1-[x(i)*x(i)];
15
end
16 17
for k=1:maxk+1 time(k) = (k-1)*dt;
18 19
end
20 21
% metode eksplisit
22
for k=2:maxk+1
23
u(1,k)=1;
24
for i=2:n+1 u(i,k) =b * u(i-1,k-1) * u(i-1,k-1) + u(i,k-1) * (1-b*u(i,k-1));
25 26
end
27 28
end
29 30
hold on;
31
plot(x,u(:,maxk/2+1),'r--');
32
plot(x,u(:,1),'blue');
33
title('explicit method')
34
xlabel('X')
35
ylabel('u(x,t)')
36
hold off
Universitas Sumatera Utara
52 Lampiran 3: Aliran Dua Fase
1
% input parameter
2
L = 1;
3
tfinal = a;
4
maxk = t;
%jumlah perputaran waktu
5
n = d;
%jumlah diskritisasi
6
dt = tfinal/maxk;
7
dx = L/n;
8
b=dt/dx;
%panjang domain
9 10
% kondisi awal
11
x(1) = 0.0;
12
u=zeros(n+1,1);
13
u(1,1) = 1.0;
14 15
for i = 2:n+1
16
x(i) =x(i-1) + dx;
17
u(i,1) = 0.0;
18
end
19 20
for k=1:maxk+1 time(k) = (k-1)*dt;
21 22
end
23 24
% metode eksplisit
25
F = zeros(n+1,1);
26
for k = 2:maxk+1
27
u(1,k)=1;
28
v = velocity(u(:,k-1),dx);
29 30
F = Fflux(u(:,k-1))*v; Fkanan = F(2:n+1,1);
31
Fkiri
u(2:n+1,k) = u(2:n+1,k-1) - b*(Fkanan - Fkiri);
32 33
= F(1:n,1);
end
34 35
% representasi grafik
36
hold on;
37
plot(x,u(:,maxk/2+1),'red');
38
plot(x,u(:,maxk+1),'black');
39
plot(x,u(:,1),'blue');
Universitas Sumatera Utara
53
40
title('Profile Saturation')
41
xlabel('x')
42
ylabel('S(x,t)')
43
hold off;
fungsi yang diimplementasikan pada code Aliran dua fase diatas
1
% subfungsi fungsi flux
2
function value = Fflux(u)
3
vr = k;
% nilai viskositas
4
value=(vr*u.*u)./(vr*u.*u+(1-u).*(1-u));
1
% subfungsi velocity
2
function vel = velocity(u,dx)
3
m = size(u);
4
laval = 0.0;
5
for j=2:m increment = (1/lamda(u(j-1))+1/lamda(u(j)))*0.5*dx; laval = increment+ laval;
6 7 8
end
9
vel = 1/laval;
fungsi lamda diimplementasikan pada subfungsi velocity
1
function val = lamda(u)
2
vr = k;
% nilai viskositas
3
val=vr*u.*u+(1-u).*(1-u);
Universitas Sumatera Utara
54 Lampiran 4: Aliran Dua Fase Nonekuilibrium
1
L = 1.0;
% Panjang domain
2
M = 400;
% jumlah diskritisasi
3
dx = L/M;
4
x = linspace(0,L,M+1)';
5 6
dt = 0.001
7 8
tau = 0.02;
9
b = dt/dx;
10
c = tau/dt;
11
cc = tau/dx;
12
maxt = 100;
% jumlah waktu
13 14
sigma = initialsigma(cc, M); % kondisi awal untuk sigma
15
17
U = zeros(2*M,1); U(1:M,1) = 0.0;
18
U(M+1:2*M,1) = sigma(2:M+1,1);
16
19 20
for k = 1:maxt+1
21
Uold = U;
22
U = iterasiU(U, Uold, b, c, dx);
23
end
24 25
sat = [1; U(1:M,1)];
26 27
effsat = [1; U(M+1:2*M,1)]; figure(1);
28
hold on;
29
plot(x,sigma,'b--');
30
plot(x,sat, 'k', x,effsat,'r');
31
legend('sigma awal', 'saturation', 'sigma');
Terdapat beberapa subfungsi yang diimplementasikan pada M-File diatas yaitu
1
%subfungsi initialsigma akan menghasilkan nilai awal sigma
2
function sigma = initialsigma(c, M)
3
Niter = 100;
Universitas Sumatera Utara
55
4
sigma = zeros(M+1,1);
5
sigma(1) = 1.0;
6
si = 0.2; %Initial Guess
7
% Implementation of Newton's Iteration
8
for i = 2:M+1
9
ss = sigma(i-1);
10
fluxfunold = fluxfunc(ss);
11
for k=1:Niter
12
fluxfun = fluxfunc(si);
13
derfluxfun = derfluxfunc(si);
14 15
Fsi=si+c*(fluxfun-fluxfunold); Fsi2=1+c*derfluxfun;
16
Delta = -Fsi/Fsi2;
17
si = si + Delta;
18
if (Delta<1.0e-8) break;
19
end
20 21
end
22
sigma(i) = si;
23
end
24
end
1
%subfungsi iterasi newton pada U
2
function [U] = iterasiU(U, Uold, b, c, dx)
3
Niter = 100;
% jumlah iterasi
4 5
for m = 1:Niter
6
F = Fvektor(U, Uold, b, c, dx);
7
jacobi = matriksjacobian(U, b, c);
8
Delta = -jacobi \ F;
9
U = U + Delta; if (norm(Delta)<1.0e-8)
10 11
m;
12
break; end
13 14
end
15
end
subfungsi iterasiU terdapat subfungsi Fvektor dan matriksjacobian berikut
Universitas Sumatera Utara
56
1
% subfungsi untuk menghitung F(U)
2
function [F] = Fvektor(U, Uold, b, c, dx)
3 4
tt = length(U);
5
F = zeros(tt,1);
6
M = tt/2;
7
v = velocity(U(1:M+1), dx);
8
F(1) = U(1) - Uold(1) + v*b*(fluxfunc(0.5*(U(1+M) + Uold(M+1)))) - (v*b*(0.5*(1+1)));
9 10
F(M+1) = 0.5*(U(M+1) - U(1) + Uold(M+1) - Uold(1)) c*(U(1)-Uold(1));
11 12 13 14
for i=2:M F(i) = U(i) - Uold(i) + v*b*(fluxfunc(0.5*(U(i+M)+ Uold(M+i)))) - v*b*(fluxfunc(0.5*(U(i+M-1)+
15 16
Uold(M+i-1))));
17
F(M+i) = 0.5*(U(M+i) - U(i) + Uold(M+i) - Uold(i)) c*(U(i)-Uold(i));
18 19 20
end
21
end
22
%
Pada subfungsi iterasiU dihitung matriks jacobian F(U) yang diimplementasikan oleh subfungsi berikut
1
%fungsi matriksjacobi
2
function [jacobian] = matriksjacobian(U, b, c)
3 4
tt = length(U);
5
M = tt/2;
6
jacobian = zeros(tt);
7
for l=1:M
8
jacobian(l,l) = 1;
9
jacobian(l+M,l+M) = 0.5;
10 11
jacobian(l+M, l) = -c-0.5; end
12 13
for l=1:M
Universitas Sumatera Utara
57
14 15
jacobian(l,l+M) = 0.5*b*derfluxfunc(U(l+M)); end
16 17
for l=2:M
18
jacobian(l,M+l-1) = -0.5*b*derfluxfunc(U(M+l-1)); end
19 20 21 22
jacobian = sparse(jacobian); end
Pada subfungsi Fvektor terdapat fungsi fluxfunc dan velocity
1
% subfungsi fungsi flux
2
function new = fluxfunc(sigma)
3
vr = 7;
4
new=(vr.*sigma.*sigma)/(vr.*sigma.*sigma+(1-sigma).*(1-sigma));
1
% subfungsi velocity sebagai nilai persamaan Darcy
2
function vel = velocity(u,dx)
3
m=size(u);
4
laval = 0.0;
5
for j=2:m increment = (1/lamda(u(j-1))+1/lamda(u(j)))*0.5*dx; laval = increment+ laval;
6 7 8 9 10
end vel=1/laval; end
pada fungsi matriksjacobi terdapat subfungsi derflux
1
%subfungsi untuk menghitung turunan fungsi flux
2
function val=derfluxfunc(sigma)
3
vr = 7;
4
val = 2*vr*sigma.*(1-sigma)/((1-sigma).ˆ2+vr*sigma.ˆ2)ˆ2;
Universitas Sumatera Utara
