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

1 vue (au cours des 30 derniers jours)
MINATI
MINATI le 29 Avr 2019
Modifié(e) : MINATI le 30 Avr 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 commentaires
David Wilson
David Wilson le 30 Avr 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

Connectez-vous pour commenter.

Réponse acceptée

David Wilson
David Wilson le 30 Avr 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 commentaire
MINATI
MINATI le 30 Avr 2019
Modifié(e) : MINATI le 30 Avr 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

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Logical dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by