Dear All
I wrote this small code to calculate unsteady PDE but I got error message. I could not solve it because I am still beginer in using matlab
The error message is:
??? Subscripted assignment dimension mismatch.
Error in ==> FDDD at 55
S(j,i)=S(j,i-1)+dt*(K1*S_D2+K2*D1_2+K33.*D1-S_1);
the Code:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Script to solve the counter-current imbibition in reservoir
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
phi=0.3;
mun=0.01;
k=0.02*10^-12;
rohn=660;
rohi=1000;
droh=rohi-rohn;
gamans=0.03;
gamais=.075;
sm=100*10^5;
beta=(gamans-gamais)*sm;
c=0.1;
g=9.8;
K1=beta*k/(phi*mun);
K2=-K1;
K3=k*g*c*droh/(phi*mun);
K4=2*K3;
% K3=0;
% K4=0;
%%% Initialize variables
dz = .01; %each depth step
Nz = 100; % Choose the number of depth steps
Nt = 5000; % Choose the number of time steps
dt =(24*60*60)/Nt; %Length of each time step in seconds
S0=0.;
S = S0*ones(Nz+1,Nt+1); %Create saturation matrix with Nz+1 rows, and Nt+1 columns
time = [0:12/Nt:12];
S(1,:) = 1; %Set surface saturation
% X axis
z0=0;
z(1) = z0;
for i = 1:Nz
z(i+1) = z0 + i*dz;
end
% z';
% z(100)';
maxiter = 500;
for iter = 1:maxiter
Slast = S; %Save the last guess
S(:,1) = Slast(:,end); %Initialize S at t=0 to the last temp
for i=2:Nt+1,
for j=2:Nz
D2 = (S(j+1,i-1)-2*S(j,i-1)+S(j-1,i-1))/dz^2;
S_D2=(1-S(j,i-1)).*D2;
D1 = (S(j,i-1)-S(j-1,i-1))/dz;
D1_2 =D1.*D1;
K33=K3./z.^2;
K44=K4./z.^3;
S_1=K44.*(S(j,i-1)-1);
S(j,i)=S(j,i-1)+dt*(K1*S_D2+K2*D1_2+K33.*D1-S_1);
S(Nz+1,i) = 0; % Enforce bottom BC
% S(Nz+1,i) = 0; % Enforce bottom BC
end
end
err(iter) = max(abs(S(:)-Slast(:))); %Find difference between last two solutions
if err(iter)<1E-4
break; % Stop if solutions very similar, we have convergence
end
end
if iter==maxiter;
warning('Convergence not reached')
end
figure(1)
plot(log(err)), title('Convergence plot')
% figure(2)
% imagesc([0 12],[0 1],S); title('Saturation plot (imagesc)')
% colorbar
figure(3)
depth = [0:dz:Nz*dz];
contourf(time,-depth,S); title('Saturation plot (contourf)')
colorbar
figure(4)
plot(time,S(1,:),time,S(21,:),time,S(41,:),time,S(61,:))
xlabel('Time (months)'); ylabel('Saturation');
legend('0m','5m','10m','15m', '20m')
|
Search |
Post Question |
Groups |
Stream |
About |
Register
Mohamed
197355 articles. 
486 followers. 
4 Replies
386 Views
