Please help me to plot surface figure. I want to draw this attached photo 2 dimension of (NG,etta) in 3 dimension as (NG, R,etta)

3 views (last 30 days)
Tarek
Tarek on 11 Jul 2025
Edited: Torsten on 11 Jul 2025
%Error
% Attempting to evaluate the solution outside the interval [0.000000e+00, 1.000000e+00] where it is defined.
%Error in proj (line 65)
% Z(i, :) = deval(sol,N,1); % Store the z-axis data
%code
function sol= proj
clc;clf;clear;
global n;
global s
%Relation of base fluid
rhof=997.1*10^-3;kf=0.613*10^5;cpf=4179*10^4;muf=10^-3*10;
alfaf=kf/(rhof*cpf);
bef=21*10^-5;
ky=muf/rhof;
disp('ky');disp((muf/rhof));
%sigf=0.05*10^-8;
%Ag
ph1=0.01;
rho1=10500*10^-3;
cp1=235*10^4;
k1=429*10^5;be1=21*10^-5;
%sig1=0.74*10^-2;
%copper
ph2=0.01;
rho2=8933*10^-3;
cp2=385*10^4;
k2=400*10^5;
%sig2=5.96*10^-1;
be2=1.67*10^-5;
%Alumina
ph3=0.01;
rho3=3970*10^-3;
cp3=765*10^4;
k3=40*10^5;
be3=0.85*10^-5;
%sig3=3.5*10^-1;
%Relation of ternary hyprid
kn=kf*((k3+2*kf-2*ph3*(kf-k3))/(k3+2*kf+ph3*(kf-k3)));
kh=kn*((k2+2*kn-2*ph2*(kn-k2))/(k2+2*kn+ph2*(kn-k2)));
kt=kh*((k1+2*kh-2*ph1*(kh-k1))/(k1+2*kh+ph1*(kh-k1)));
mut= muf/((1-ph1)^2.5*(1-ph2)^2.5*(1-ph3)^2.5);
rhot=(1-ph1)*((1-ph2)*((1-ph3)+ph3*(rho3/rhof))+ph2*(rho2/rhof))+ph1*(rho1/rhof);
%vt=rhot*cpt
vt =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*cp3)/(rhof*cpf)))+ph2*((rho2*cp2)/(rhof*cpf)))+ph1*((rho1*cp1)/(rhof*cpf));
%disp('vt');disp(vt);
%vb=rho*betb
vb =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*be3)/(rhof*bef)))+ph2*((rho2*be2)/(rhof*bef)))+ph1*((rho1*be1)/(rhof*bef));
%disp('vb');disp(vb);disp(ky);
myLegend1 = {};myLegend2 = {};
rr = [1 2 3 4]
numR = numel(rr);
m = linspace(0,1);
a=0.001;b=0.001;p=.001/((1-0.01)*(mut/muf)*(rhof/rhot));
Ec=10;R=0.55;at=0.01;gamma=pi/4;
prf=6.9;d=0.009;e=0.1;Rd=0.45;
Tw=273+50;Ti=273+27;deltaT=Tw-Ti;
disp('coe');disp((mut/muf)*(rhof/rhot));
Lf=rhof*kf;
y0 = [1,0,1,0,0,1,0,1];options =bvpset('stats','on','RelTol',1e-5);
Z = zeros(numR, length(m));
for i = 1:numR
R= rr(i);
solinit = bvpinit(m, y0);
sol = bvp4c(@projfun, @projbc, solinit, options);
N=(kt/kf)*(1+Rd)*d*R*(1/(muf/rhof)^2)*(sol.y(7,:).^2)+((mut/kf)*(1/Ti)*(1/at^2))*R*(d*(muf/rhof)*(sol.y(2,:).^2)+e*(sol.y(4,:).^2)+R*d*p*(sol.y(1,:).^2));
Z(i, :) = deval(sol,N,1); % Store the z-axis data
end
[X, Y] = meshgrid(m, rr);
shading interp;
surf(X, Y, N);
xlabel('eta');
ylabel('Gr');
zlabel('NG');
title('Variation of velocity with Grashof number,Gr in 3D' );
grid on
shading flat;
colorbar;
function dy= projfun(~,y)
dy= zeros(8,1);
% alignComments
E = y(1);
dE = y(2);
F = y(3);
dF= y(4);
w = y(5);
dw=y(6);
t = y(7);
dt = y(8);
dy(1) = dE;
dy(2) = (((rhot/mut)*(a*(muf/rhof)^0.5*(E*F+E^2)+a*(muf/rhof)*w*dE-(mut/muf)*(rhof/rhot)*p*(1-0.01)*E+R*a*(muf/rhof)*sin(gamma)*(vb/(rhof*bef))*t)));
dy(3) = dF;
dy(4) = (((rhot/mut)*(b*(muf/rhof)^0.5*(F^2+F*E)+(muf/rhof)*b^0.5*a^(1.5)*dF)));
dy(5) =-(a*F+b*E);
dy(6) = (((rhot/mut)*((muf/rhof)^0.5*w*dw+R*b*(muf/rhof)*cos(gamma)*(vb/(rhof*bef))*t)));
dy(7) = dt;
dy(8)=prf*(1/(kt/kf))*(1/(1+((prf*Rd)/((kt/kf)))))*((vt/(rhof*cpf))*(muf/rhof)^0.5*w*dt-(mut/muf)*Ec*1*dw^2) ;
end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
ya(3)-1;
ya(5);
ya(6);
ya(7)+1-(1/0.9)*ya(8);
yb(1)-0.01;
yb(3);
yb(7);
% yb(7);
];
end

Sign in to answer this question.

Accepted Answer

Torsten
Torsten on 11 Jul 2025
Edited: Torsten on 11 Jul 2025
Ran in:
proj()
ky 0.0100
rr = 1×4
1 2 3 4
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
coe 0.8935 The solution was obtained on a mesh of 285 points. The maximum residual is 2.085e-06. There were 15848 calls to the ODE function. There were 185 calls to the BC function. The solution was obtained on a mesh of 308 points. The maximum residual is 1.901e-06. There were 15138 calls to the ODE function. There were 105 calls to the BC function. The solution was obtained on a mesh of 352 points. The maximum residual is 4.438e-06. There were 13545 calls to the ODE function. There were 80 calls to the BC function. The solution was obtained on a mesh of 352 points. The maximum residual is 7.622e-06. There were 9852 calls to the ODE function. There were 54 calls to the BC function.
ans = struct with fields:
solver: 'bvp4c' x: [0 0.0076 0.0152 0.0227 0.0303 0.0370 0.0438 0.0513 0.0589 0.0665 0.0741 0.0816 0.0892 0.0968 0.1044 0.1120 0.1195 0.1271 0.1347 0.1423 0.1498 0.1574 … ] (1×352 double) y: [8×352 double] yp: [8×352 double] stats: [1×1 struct]
function sol= proj
clc;clf;clear;
global n;
global s
%Relation of base fluid
rhof=997.1*10^-3;kf=0.613*10^5;cpf=4179*10^4;muf=10^-3*10;
alfaf=kf/(rhof*cpf);
bef=21*10^-5;
ky=muf/rhof;
disp('ky');disp((muf/rhof));
%sigf=0.05*10^-8;
%Ag
ph1=0.01;
rho1=10500*10^-3;
cp1=235*10^4;
k1=429*10^5;be1=21*10^-5;
%sig1=0.74*10^-2;
%copper
ph2=0.01;
rho2=8933*10^-3;
cp2=385*10^4;
k2=400*10^5;
%sig2=5.96*10^-1;
be2=1.67*10^-5;
%Alumina
ph3=0.01;
rho3=3970*10^-3;
cp3=765*10^4;
k3=40*10^5;
be3=0.85*10^-5;
%sig3=3.5*10^-1;
%Relation of ternary hyprid
kn=kf*((k3+2*kf-2*ph3*(kf-k3))/(k3+2*kf+ph3*(kf-k3)));
kh=kn*((k2+2*kn-2*ph2*(kn-k2))/(k2+2*kn+ph2*(kn-k2)));
kt=kh*((k1+2*kh-2*ph1*(kh-k1))/(k1+2*kh+ph1*(kh-k1)));
mut= muf/((1-ph1)^2.5*(1-ph2)^2.5*(1-ph3)^2.5);
rhot=(1-ph1)*((1-ph2)*((1-ph3)+ph3*(rho3/rhof))+ph2*(rho2/rhof))+ph1*(rho1/rhof);
%vt=rhot*cpt
vt =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*cp3)/(rhof*cpf)))+ph2*((rho2*cp2)/(rhof*cpf)))+ph1*((rho1*cp1)/(rhof*cpf));
%disp('vt');disp(vt);
%vb=rho*betb
vb =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*be3)/(rhof*bef)))+ph2*((rho2*be2)/(rhof*bef)))+ph1*((rho1*be1)/(rhof*bef));
%disp('vb');disp(vb);disp(ky);
myLegend1 = {};myLegend2 = {};
rr = [1 2 3 4]
numR = numel(rr);
m = linspace(0,1);
a=0.001;b=0.001;p=.001/((1-0.01)*(mut/muf)*(rhof/rhot));
Ec=10;R=0.55;at=0.01;gamma=pi/4;
prf=6.9;d=0.009;e=0.1;Rd=0.45;
Tw=273+50;Ti=273+27;deltaT=Tw-Ti;
disp('coe');disp((mut/muf)*(rhof/rhot));
Lf=rhof*kf;
y0 = [1,0,1,0,0,1,0,1];options =bvpset('stats','on','RelTol',1e-5);
Z = zeros(numR, length(m));
for i = 1:numR
R= rr(i);
if i==1
solinit = bvpinit(m, y0);
else
guess = @(x)interp1(sol.x,(sol.y).',x);
solinit = bvpinit(sol.x,guess);
end
sol = bvp4c(@projfun, @projbc, solinit, options);
N=(kt/kf)*(1+Rd)*d*R*(1/(muf/rhof)^2)*(sol.y(7,:).^2)+((mut/kf)*(1/Ti)*(1/at^2))*R*(d*(muf/rhof)*(sol.y(2,:).^2)+e*(sol.y(4,:).^2)+R*d*p*(sol.y(1,:).^2));
Z(i,:) = interp1(sol.x,N,m);
end
[X, Y] = meshgrid(m, rr);
shading interp;
surf(X, Y, Z);
xlabel('eta');
ylabel('Gr');
zlabel('NG');
title('Variation of velocity with Grashof number,Gr in 3D' );
grid on
shading flat;
colorbar;
function dy= projfun(~,y)
dy= zeros(8,1);
% alignComments
E = y(1);
dE = y(2);
F = y(3);
dF= y(4);
w = y(5);
dw=y(6);
t = y(7);
dt = y(8);
dy(1) = dE;
dy(2) = (((rhot/mut)*(a*(muf/rhof)^0.5*(E*F+E^2)+a*(muf/rhof)*w*dE-(mut/muf)*(rhof/rhot)*p*(1-0.01)*E+R*a*(muf/rhof)*sin(gamma)*(vb/(rhof*bef))*t)));
dy(3) = dF;
dy(4) = (((rhot/mut)*(b*(muf/rhof)^0.5*(F^2+F*E)+(muf/rhof)*b^0.5*a^(1.5)*dF)));
dy(5) =-(a*F+b*E);
dy(6) = (((rhot/mut)*((muf/rhof)^0.5*w*dw+R*b*(muf/rhof)*cos(gamma)*(vb/(rhof*bef))*t)));
dy(7) = dt;
dy(8)=prf*(1/(kt/kf))*(1/(1+((prf*Rd)/((kt/kf)))))*((vt/(rhof*cpf))*(muf/rhof)^0.5*w*dt-(mut/muf)*Ec*1*dw^2) ;
end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
ya(3)-1;
ya(5);
ya(6);
ya(7)+1-(1/0.9)*ya(8);
yb(1)-0.01;
yb(3);
yb(7);
% yb(7);
];
end

More Answers (0)

Categories

Find more on Gravitation, Cosmology & Astrophysics in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!