Dear experts,
Anyone help me,
how to plot isotherm lines, contour plot, surface plot using matlab bvp4c.
suggests any sample codes or any ideas.
Thank you.
this is our code
function hybrid(solver)
if nargin < 1
solver = 'bvp4c';
end
bvpsolver = fcnchk(solver);
xlow=0;
g1=-1;
g2=1;
g3=0.5;
g4=0.5;
global St rad m
St=0.2;
infinity=01;
maxinfinity=15;
solinit=bvpinit(linspace(0,infinity,1001),[St g1 1 g2 0 1 g3 1 g4]);
sol = bvpsolver(@ode,@bc,solinit);
x = sol.x;
y = sol.y;
fprintf('\n');
fprintf('At eta= %g \t %10.8f \t %10.8f \t %10.8f %10.8f.\n',xlow,y(2,1),y(4,1),y(7,1),y(9,1))
for Bnew = infinity:01:maxinfinity
solinit = bvpxtend(sol,Bnew);
sol = bvpsolver(@ode,@bc,solinit);
x = sol.x;
y = sol.y;
fprintf('At eta= %g \t %10.8f \t %10.8f \t %10.8f %10.8f.\n',Bnew,y(2,1),y(4,1),y(7,1),y(9,1))
end
fprintf('%10.8f \t %10.8f \t %10.8f \t %10.8f.\n',y(2,1),y(4,1),y(7,1),y(9,1))
plot(x,y(6,:),'-','Linewidth',1.5)
% hold on
%Cf=rad*(1+rad)^m*((0.32906655)^2+(-1.12865255)^2)^(1/2)
%Nu=-rad*(1+rad)^m*y(7);
%Sh=-rad*(1+rad)^m*y(9)
%Sn=-rad*(1+rad)^m*y(11)
function dydx=ode(x,y)
% rhos1 = 2810; Cps1 = 960; ks1 = 173; Pi1 = 0.1; ss1 = 26.77*(10^6); %AA7075
% rhos2 = 2720; Cps2 = 893.2; ks2 = 222; Pi2 = 0.3; ss2 = 34.83*(10^6); %AA7072
% rhof = 792; Cpf = 2545; kf = 0.2035; sf = 0.5*(10^(-6)); %(Methonal)
Np=2;
alpha=0.15;
eps=0.45;
Re=1;
Pr=1.5;
rad=0.2;
m=1.1;
n=0.1;
Sc=01;
Beta=1;
B=1;
Ha=0.2;
Ec=0.1;
Del=0.01;
h=0.5;
Rp=1.5;
% knf = kf*(ks1+(Np-1)*kf-(Np-1)*Pi1*(kf-ks1))/(ks1+(Np-1)*kf+Pi1*(kf-ks1));
% K =(ks2+(Np-1)*knf+Pi2*(knf-ks2)) /(ks2+(Np-1)*knf-(Np-1)*Pi2*(knf-ks2));%khnf/kf
%
% snf = sf*(ss1+2*sf-2*Pi1*(sf-ss1))/(ss1+2*sf+Pi1*(sf-ss1));
% S = (ss2+2*snf-2*Pi2*(snf-ss2))/(ss2+2*snf+Pi2*(snf-ss2));%shnf/sf
%
% Rh=(1-Pi2)*((1-Pi1)+Pi1*((rhos1*Cps1)/(rhof*Cpf)))+Pi2*((rhos2*Cps2)/(rhof*Cpf));
% Mu=1/((1-Pi1)^(2.5)*(1-Pi2)^(2.5));
% D=1.8;
% RaggRp=3.34;
% Pmax=0.605;
% Pint=(RaggRp)^(D-3);
% Pagg=Pi/Pint;
% rhoagg=(1-Pint)*rhof+Pint*rhos;
% Mu=(1/(Pagg/Pmax))^(-2.5*Pmax);
% Rh=(1-Pagg)+((Pagg*rhoagg)/rhof);
% rhoCpagg=(1-Pint)*(rhof*Cpf)+Pint*rhos*Cps;
% kagg=(kf/4)*((3*Pint-1)*(ks/kf)+(3*(1-Pint)-1)+(((3*Pint-1)*(ks/kf)+(3*(1-Pint)-1))^2+(8*ks/kf))^0.5);
% RhCp=(1-Pagg)+Pagg*(rhoCpagg/(rhof*Cpf));
% K=(kagg+2*kf+2*Pagg*(kagg-kf))/(kagg+2*kf-Pagg*(kagg-kf));
% S=1+(((3*(ss/sf)-1)*Pagg)/((ss/sf)+2)*((ss/sf)-1)*Pagg);
rhos = 2810; Cps = 960; ks = 173; Pi = 0.1; ss = 26.77*(10^6); %AA7075
rhof = 792; Cpf = 2545; kf = 0.2035; sf = 0.5*(10^(-6)); %(Methonal)
Mu=1+2.5*Del+4.5*(1/((h/Rp)*(2+(h/Rp))*(1+(h/Rp))^2));
Rh=(1-Del)+Del*(rhos/rhof);
RhCp=(1-Del)+Del*(rhos*Cps/rhof*Cpf);
K=(ks+2*kf-2*Del*(kf-ks))/(ks+2*kf+Del*(kf-ks));
S=(ss+2*sf-2*Del*(sf-ss))/(ss+2*sf+Del*(sf-ss));
dydx=[y(2)
Rh*(y(1)^2-y(3)^2+m*eps*(x+alpha)*y(1)*y(2)+y(5)*y(2))-(B/Re^2)+S*Ha*y(1)/Mu*Re^((1-n)/(1+n))*(1+rad)^(2*m)*(1+(1/Beta))
y(4)
Rh*(2*y(1)*y(3)+y(5)*y(4)+m*eps*(x+alpha)*y(1)*y(4))+S*Ha*y(3)/Mu*Re^((1-n)/(1+n))*(1+rad)^(2*m)*(1+(1/Beta))
-2*y(1)-m*eps*(x+alpha)*y(2)
y(7)
(y(1)*y(7)*m*eps*(x+alpha)*Pr+Pr*y(5)*y(7)+Pr*(S*RhCp)*Ha*Ec*(y(1)^2+y(3)^2))/(K*RhCp)*Re^((1-n)/(1+n))*(1+rad)^(2*m)
y(9)
(Sc*y(1)*y(9)*m*eps*(x+alpha)+Sc*y(5)*y(9))/(Mu*Re^((1-n)/(1+n))*(1+rad)^(2*m))
];
end
function res=bc(y0,yinf)
res=[y0(1)-St
y0(3)-1
y0(5)
y0(6)-1
y0(8)-1
yinf(1)
yinf(3)
yinf(6)
yinf(8)
];
end
end