No error but I am getting Complex double!

12 次查看(过去 30 天)
clear ALL
close ALL
clc
Area=0.1963;
L=4;
porosity=0.18;
K=61.5;
injrate=0.01;
oilvisc=4.2;
watervisc=0.85;
Boi=1.21;
Swr=0.28;
Sorw=0.35;
KrwMax=0.31;
KroMax=0.65;
Nw=2.25;
No=2.45;
CO=0.000001;
CW=0.000004;
CR=0.000003;
POi=1500;
PWF=1460;
dt=0.2;
Vr=Area*0.4;
for i=1:10
Swi(1,i)=Swr;
PO(1,i)=POi;
PCOW(1,i)=5;
Krw(1,i)=KrwMax*((Swi(1,i)-Swr)/(1-Sorw-Swr))^Nw;
Kro(1,i)=KroMax*((1-Swi(1,i)-Sorw)/(1-Sorw-Swr))^No;
Omob(1,i)= Kro(1,i)/oilvisc;
Wmob(1,i)=Krw(1,i)/watervisc;
end
for n=1:15
for i=1:10
Ctot(n,:)=CR+Swi(n,:)*CW+(1-Swi(n,:))*CO;
if i==1
Ttot=(2.6366e-4)*K*(Wmob(n,i)+Omob(n,i))*(Area/0.4);
Tw=(2.6366e-4)*K*(Wmob(n,i))*(Area/0.4);
A(i,i)=-(Ttot+(Vr*porosity*Ctot(n,i))/(dt)); %E
A(i,i+1)=Ttot; %F
R(i,1)= (-(Vr*porosity*Ctot(n,i)*PO(n,i))/(dt))+Tw*(PCOW(n,i+1)-PCOW(n,i))-injrate; %R
elseif i==10
WI(n,i)=(2.6366e-4)*K*(Wmob(n,i)+Omob(n,i))*(Area/(0.4/2));
TtotL=Ttot;
TwL=Tw;
Ttot=0;
Tw=0;
%Ttot=2.6366*10^-4*K*(Wmob(n,i)+Omob(n,i))*(A/0.4);
%Tw=2.6366*10^-4*K*(Wmob(n,i))*(A/0.4);
A(i,i-1)=TtotL; %D
A(i,i)=-(TtotL+((Vr*porosity*Ctot(n,i))/(dt))+WI(n,i)); %E
R(i,1)=(-(Vr*porosity*Ctot(n,i)*PO(n,i))/(dt))-TwL*(PCOW(n,i)-PCOW(n,i-1))-WI(n,i)*PWF; %R
else
TtotL=Ttot;
TwL=Tw;
Ttot=(2.6366e-4)*K*(Wmob(n,i)+Omob(n,i))*(Area/0.4);
Tw=(2.6366e-4)*K*(Wmob(n,i))*(Area/0.4);
A(i,i)=-(TtotL+Ttot+(Vr*porosity*Ctot(n,i))/dt); %E
A(i,i+1)=Ttot; %F
A(i,i-1)=TtotL; %D
R(i,1)= (-(Vr*porosity*Ctot(n,i)*PO(n,i))/(dt))+Tw*(PCOW(n,i+1)-PCOW(n,i))-TwL*(PCOW(n,i)-PCOW(n,i-1)); %R
end
end
PO(n+1,:)=A\R;
for i=1:10
if i==1
q_t(n,i)=injrate;
Tw=(2.6366e-4)*K*(Wmob(n,i))*(Area/0.4);
Swi(n+1,i)= (1+(CR+CW)*(PO(n+1,i)-PO(n,i)))*Swi(n,i)+...
(dt/(Vr*porosity))*(Tw*(PO(n+1,i+1)-PO(n+1,i))-...
Tw*(PCOW(n,i+1)-PCOW(n,i))+q_t(n,i));
elseif i==10
TwL=Tw;
Tw=0;
WI(n,i)=(2.6366e-4)*K*(Wmob(n,i)+Omob(n,i))*(Area/(0.4/2));
qt(n,i)= -WI(n,i)*(PO(n+1,i)-PWF);
Swi(n+1,i)=(1+(CR+CW)*(PO(n+1,i)-PO(n,i)))*Swi(n,i)+dt/(Vr*porosity)*...
(-TwL*(PO(n+1,i)-PO(n+1,i-1))+TwL*(PCOW(n,i)-PCOW(n,i-1))+qt(n,i));
else
TwL=Tw;
Tw=2.6366e-4*K*(Wmob(n,i))*(Area/0.4);
Swi(n+1,i)=(1+(CR+CW)*(PO(n+1,i)-PO(n,i)))*Swi(n,i)+ ...
(0.2/(Vr*porosity))*Tw*(PO(n+1,i+1)-PO(n+1,i))-...
TwL*(PO(n+1,i)-PO(n+1,i-1))-Tw*(PCOW(n,i+1)-PCOW(n,i))+TwL*(PCOW(n,i)-PCOW(n,i-1));
end
if Swi(n+1,i)<Swr
Swi(n+1,i)=Swr;
end
end
for i=1:10
if Swi(n+1,i)<=Swr
PCOW(n+1,i)=5;
elseif Swi(n+1,i)>=(1-Sorw)
PCOW(n+1,i)=-5;
elseif Swi(n+1,i)>Swr && Swi(n+1,i)<0.5
PCOW(n+1,i)= min(-0.8*log((Swi(n+1,i)-Swr)/(0.5-Swr)),5);
elseif Swi(n+1,i)>0.5 && Swi(n+1)<(1-Sorw)
PCOW(n+1,i)=max(0.8*log((1-Swi(n+1,i)-Sorw)/(1-0.5-Sorw)),-5);
end
end
for i=1:10
Krw(n+1,i)=KrwMax*((Swi(n+1,i)-Swr)/(1-Sorw-Swr))^Nw;
Kro(n+1,i)=KroMax*((1-Swi(n+1,i)-Sorw)/(1-Sorw-Swr))^No;
Omob(n+1,i)= Krw(n+1,i)/watervisc;
Wmob(n+1,i)=Kro(n+1,i)/oilvisc;
end
end

回答(2 个)

Mark Sherstan
Mark Sherstan 2019-3-1
You should not use i or j as the names of loop variables as these are both names for the inbuilt imaginary unit. Try changing that and see if the problem persisits.
  1 个评论
Sarah A Alruwayi
Sarah A Alruwayi 2019-3-1
My friend has exactly the same code and got it right and she named the loops i as well.So, I do not think this is the problem

请先登录,再进行评论。


Walter Roberson
Walter Roberson 2019-3-1
When n = 4, i = 1, then 1-Swi(n+1),i)-Sorw becomes negative, so (1-Swi(n+1,i)-Sorw)/(1-Sorw-Swr) becomes negative. You are raising a negative value to a non-integer. A^B with scalar A, B is defined as exp(B*LOG(A)) so when A<0 the log is complex and unless B is an integer you get a complex result from the ^ operation. Everything else gets tainted from there.

类别

Help CenterFile Exchange 中查找有关 Logical 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by