MATLAB Answers

MINATI
0

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

Asked by MINATI
on 29 Apr 2019
Latest activity Edited by 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
[EDITED, Jan, Attachment added].

  11 Comments

@Walter
Relevant equations are in (7) and boundary conditions are in (8).
Yes we need numeric solution.
But can it be possible to derive theoretical solution from the code.
whether it is a "stiff" system. I dont know what is "stiff".
does that hint that we could provide the Jacobian or Hessian as functions?
Ans: I think its NO
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.

1 Answer

Answer by David Wilson on 30 Apr 2019
 Accepted Answer

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

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.