Asked by MINATI
on 29 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].

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.

Opportunities for recent engineering grads.

Apply Today
## 11 Comments

## Walter Roberson (view profile)

## 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

## MINATI (view profile)

## 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

## Jan (view profile)

## 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

## MINATI (view profile)

## 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

## MINATI (view profile)

## 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

## Jan (view profile)

## 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

## MINATI (view profile)

## 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

## Walter Roberson (view profile)

## 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

## MINATI (view profile)

## 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

## David Wilson (view profile)

## 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

## Jan (view profile)

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