Please help me to run the surface figure from the attached code

6 ビュー (過去 30 日間)
Tarek
Tarek 2025 年 7 月 23 日
回答済み: Star Strider 2025 年 7 月 23 日
%Array indices must be positive integers or logical values.
%Error in Tarek7_2025 (line 66)
% 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));
%code
function sol= proj
clc;clf;clear;
%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.09;b=0.001;p=.001/((1-0.01)*(mut/muf)*(rhof/rhot));
Ec=10;Gr=0.55;at=0.09;gamma=pi/4;
Rd=1;
prf=6.9;d=0.0001;e=0.001;Rd=0.8;
Tw=273+50;Ti=273+27;deltaT=Tw-Ti;
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+Gr*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+Gr*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

サインインしてこの質問に回答する。

採用された回答

Star Strider
Star Strider 2025 年 7 月 23 日
Ran in:
The error is because you are missing a multiplication operator
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)));
HERE ^
and
dy(6) = (((rhot/mut)*((muf/rhof)^0.5*w*dw+R*b*(muf/rhof)*cos(gamma)*(vb/(rhof*bef))*t)));
HERE ^
and
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));
HERE ^ ^ ^
I inserted them (there are 3 in 'N').
The result otherwise appears to be appropriate.
Running your code --
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]
%rray indices must be positive integers or logical values.
%Error in Tarek7_2025/projfun (line 92)
% 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)));
%code
function sol= proj
clc;clf;clear;
%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;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
.

その他の回答 (0 件)

カテゴリ

Help Center および File Exchange で Eigenvalue Problems についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by