# How to draw Fig. 1 from the attached pdf with this code

2 views (last 30 days)
MINATI on 29 Apr 2019
Edited: MINATI on 30 Apr 2019
function main
Pr=1; G=0.1;
% phi=input('phi='); %%0,.05, .1, .15, .2
phi=0.0;
rhof=997.1;Cpf=4179;kf=0.613; %for WATER
rhos=6320;Cps=531.8;ks=76.5; %for CuO
a1=((1-phi)^2.5)*(1-phi+phi*(rhos/rhof));
a2=(1-phi+phi*((rhos*Cps)/(rhof*Cpf)));
A=(ks+2*kf+phi*(kf-ks))/(ks+2*kf-2*phi*(kf-ks)); %%%%Knf
xa=0;xb=6;
solinit=bvpinit(linspace(xa,xb,101),[0 1 0 1 0]);
sol=bvp4c(@ode,@bc,solinit);
xint=linspace(xa,xb,101);
sxint=deval(sol,xint);
figure(1)
plot(xint,(1-phi)^-2.5*sxint(3,:),'-','Linewidth',1.5); %for f''(0)/(1-phi)^2.5 vs phi
xlabel('\eta');
ylabel('f''(0)/(1-phi)^2.5');
hold on
function res=bc(ya,yb)
res=[ya(1); ya(2)-1-G*ya(3); ya(4)-1; yb(2); yb(4)];
end
function dydx=ode(x,y)
dydx=[y(2); y(3); a1*(y(2)^2-y(3)*y(1)); y(5); -A*Pr*a2*y(1)*y(5)];
end
end
Jan on 30 Apr 2019

David Wilson on 30 Apr 2019
I didn't bother draw the other 3 lines, but you just ned to make the necessary changes to gamma for that.
If you run something like what you had originally, you only want the fist point of f''().
Pr=6.2; G=0.1;
% phi=input('phi='); %%0,.05, .1, .15, .2
phi=0.0;
rhof=997.1;Cpf=4179;kf=0.613; %for WATER
rhos=6320;Cps=531.8;ks=76.5; %for CuO
a1=((1-phi)^2.5)*(1-phi+phi*(rhos/rhof));
a2=(1-phi+phi*((rhos*Cps)/(rhof*Cpf)));
A=(ks+2*kf+phi*(kf-ks))/(ks+2*kf-2*phi*(kf-ks)); %%%%Knf
BCres= @(ya,yb) ...
[ya(1); ya(2)-1-G*ya(3); ya(4)-1; yb(2); yb(4)];
fODE = @(x,y) ...
[y(2); y(3); a1*(y(2)^2-y(3)*y(1)); y(5); -A*Pr*a2*y(1)*y(5)];
xa=0;xb=8;
solinit=bvpinit(linspace(xa,xb,101),[0 1 0 1 0]);
sol=bvp4c(fODE,BCres,solinit);
xint=linspace(xa,xb,101);
sxint=deval(sol,xint);
figure(1)
plot(xint,(1-phi)^-2.5*sxint(3,:),'-','Linewidth',1.5); %for f''(0)/(1-phi)^2.5 vs phi
xlabel('\eta');
ylabel('f''(0)/(1-phi)^2.5');
Now you have to re-run the above, but change phi over the range given in the Fig.
xa=0;xb=8;
phiv = [0:0.04:0.2]';
p = []; % collect points here
for i=1:length(phiv)
phi = phiv(i);
a1=((1-phi)^2.5)*(1-phi+phi*(rhos/rhof));
a2=(1-phi+phi*((rhos*Cps)/(rhof*Cpf)));
A=(ks+2*kf+phi*(kf-ks))/(ks+2*kf-2*phi*(kf-ks)); %%%%Knf
fODE = @(x,y) ...
[y(2); y(3); a1*(y(2)^2-y(3)*y(1)); y(5); -A*Pr*a2*y(1)*y(5)];
solinit=bvpinit(linspace(xa,xb,101),[0 1 0 1 0]);
sol=bvp4c(fODE,BCres,solinit);
p(i,1) = (1-phi)^-2.5*sxint(3,1)
end
plot(phiv, p,'o-')
xlabel('\phi'); ylabel('f''''(0) & stuff')
Resultant plot is as above.
MINATI on 30 Apr 2019
Edited: MINATI on 30 Apr 2019
Many many thanks David
It worked
Where to put the loop for Gamma G=[0 0.1 0.2 0.3] which will vary the fig