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].

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.

Opportunities for recent engineering grads.

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

Start Hunting!
## 11 Comments

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699364

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699364

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699597

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699597

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699601

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699601

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699602

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699602

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699604

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699604

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699652

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699652

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699689

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699689

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699712

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699712

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699852

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699852

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699972

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699972

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699990

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#comment_699990

Sign in to comment.