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

1 view (last 30 days)
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
[EDITED, Jan, Attachment added].
  11 Comments
David Wilson
David Wilson on 30 Apr 2019
If I understand correctly, the BVP you are trying to solve has BCs at infinity. You have chosen 6 (& the paper uses 8), so you might like to validate that approximation.
What exactly is the problem ?
My solution is for the Pr=6.2, \gamma=0.1. This seems to follow your Fig 1. tmp.png

Sign in to comment.

Accepted Answer

David Wilson
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.
  1 Comment
MINATI
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

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!