I want 4 solid lines for Fr paramter . but im getting single line why?

1 次查看(过去 30 天)
function S_a
clc
clear All
%%%%note%%%%%%%%%
%%%%%%%%%%%%%%%%% decreasing behaviour with Fr=0.4,0.8,1.2,1.6,2 %%%%%%%%%%%
global Fr Pr f2 f1 F2 shi2 shi3 lambda zeta Sk u Re Kf B A H1 H2 H3 C1 C2 C3 P1 P2 P3 K1 K2 K3 Kb n Z FR
FR = 0.4:0.4:2
for i=1:numel(FR)
Fr = FR(i); %variation paramter
end
Pr = 6.2 ;
B = 0.9;
F2 = 0:0.01:0.15
for i=1:numel(F2)
f2 = F2(i);
f1 = 0.05 ;
lambda = 0.4;
zeta=0.5;
A= 0.7;
n=3.7;
Re=0.3;
C1=765;
P1=3970;
K1=40;
%%%%%%%%%%%%%%%%%%%%%%%
C2=385; % specific heat
P2=8933; % density
K2=400; % thermal conductivity
%%%%%%%%%%% %%%%%%%%%%%%
C3=4180; % specific heat
P3=997.1; % density
K3=0.6071; % thermal conductivity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%
Z=zeta+B;
H1=P1*C1; % pho*cp nanoparticle 1
H2=P2*C2; % pho*cp nanoparticle 2
H3=P3*C3; % pho*cp base fluid
Kb = (K2 + (n - 1)*K3 - f2*(n - 1)*(K3 - K2))/(K2 + (n - 1)*K3 + f2*(K3 - K2)); % khnf/kbf
Kf=((K2 + (n - 1)*K3 - f2*(n - 1)*(K3 - K2))/(K2 + (n - 1)*K3 + f2*(K3 - K2)))*((K1 + (n - 1)*K3 - f1*(n - 1)*(K3 - K1))/(K1 + (n - 1)*K3 + f1*(K3 - K1)));
% khnf/kf
shi2=(((1 - f1)^2.5)*((1 - f2)^2.5))*(((1 - f1)*(1 - f2))+f1*(P1/P3)) + f2*(P2/P3) ;
% shi2 in velocity equation
shi3=(1 - f2)*((1 - f1) + f1*((H1)/(H3))) + f2*((H2)/(H3));
u=1/((1-f1)^2.5);
options = bvpset('RelTol',1e-6,'Stats','off') ;
solinit = bvpinit(linspace(0,5,10),[1; 0 ;0 ;0 ;0; 0]);
sol = bvp4c(@odes,@bcs,solinit,options);
Sk(i)=u*1/((Re)^0.5)*(sol.y(3,1)-sol.y(2,1)/B);
end
plot(F2,Sk,'-y')
%--------------------------------------------------------------------------
function dydx = odes(~,y)
% global zeta Pr psi2 psi3 k B Fr A
dydx = zeros(6,1);
dydx(1) = y(2); % y(1)=f
dydx(2) = y(3);
dydx(3) = y(4);
dydx(4) = -(2/(zeta+B)*y(4))+(1/((zeta+B)^2)*y(3))-(1/((zeta+B)^3)*y(2))+2*lambda*(y(3)+(1/(zeta+B)*y(2)))+shi2*(-((B/(zeta+B)^2)*y(1)*y(3))-((B/(zeta+B)*y(1)*y(4)))+(B/((zeta+B)^3)*y(1)*y(2))+(3*B/((zeta+B)^2)*y(2)*y(2))+(3*B/((zeta+B)^2)*y(2)*y(3))+2*Fr*(2*y(2)*y(3)+((1/(zeta+B))*y(2)*y(2))));
dydx(5) = y(6); % y(5)=theta
dydx(6) = -(1/(zeta+B))*y(6)-((Pr*shi3*B)/(Kb*(zeta+B)))*(y(1)*y(6)-A*y(2)*y(5));% k=khnf/kbf
%--------------------------------------------------------------------------
end
function res = bcs(ya,yb)
res = [ya(1); ya(2)-1;ya(5)-1; yb(2); yb(3); yb(5)];
end
%%%%%%%%%%%%2nd%%%%%%%%%%%%%
Fr = 0.8;
Pr = 6.2 ;
B = 0.9;
F2 = 0:0.01:0.15
for i=1:numel(F2)
f2 = F2(i);
f1 = 0.05 ;
lambda = 0.4;
zeta=0.5;
A= 0.7;
n=3.7;
Re=0.3;
C1=765;
P1=3970;
K1=40;
%%%%%%%%%%%%%%%%%%%%%%%
C2=385; % specific heat
P2=8933; % density
K2=400; % thermal conductivity
%%%%%%%%%%% %%%%%%%%%%%%
C3=4180; % specific heat
P3=997.1; % density
K3=0.6071; % thermal conductivity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%
Z=zeta+B;
H1=P1*C1; % pho*cp nanoparticle 1
H2=P2*C2; % pho*cp nanoparticle 2
H3=P3*C3; % pho*cp base fluid
Kb = (K2 + (n - 1)*K3 - f2*(n - 1)*(K3 - K2))/(K2 + (n - 1)*K3 + f2*(K3 - K2)); % khnf/kbf
Kf=((K2 + (n - 1)*K3 - f2*(n - 1)*(K3 - K2))/(K2 + (n - 1)*K3 + f2*(K3 - K2)))*((K1 + (n - 1)*K3 - f1*(n - 1)*(K3 - K1))/(K1 + (n - 1)*K3 + f1*(K3 - K1)));
% khnf/kf
shi2=(((1 - f1)^2.5)*((1 - f2)^2.5))*(((1 - f1)*(1 - f2))+f1*(P1/P3)) + f2*(P2/P3) ;
% shi2 in velocity equation
shi3=(1 - f2)*((1 - f1) + f1*((H1)/(H3))) + f2*((H2)/(H3));
u=1/((1-f1)^2.5);
options1 = bvpset('RelTol',1e-6,'Stats','off') ;
solinit1 = bvpinit(linspace(0,5,10),[1; 0 ;0 ;0 ;0; 0]);
sol = bvp4c(@odes1,@bcs1,solinit1,options1);
Sk(i)=u*1/((Re)^0.5)*(sol.y(3,1)-sol.y(2,1)/B);
end
plot(F2,Sk,'--r')
%--------------------------------------------------------------------------
function dydx = odes1(~,y)
% global zeta Pr psi2 psi3 k B Fr A
dydx = zeros(6,1);
dydx(1) = y(2); % y(1)=f
dydx(2) = y(3);
dydx(3) = y(4);
dydx(4) = -(2/(zeta+B)*y(4))+(1/((zeta+B)^2)*y(3))-(1/((zeta+B)^3)*y(2))+2*lambda*(y(3)+(1/(zeta+B)*y(2)))+shi2*(-((B/(zeta+B)^2)*y(1)*y(3))-((B/(zeta+B)*y(1)*y(4)))+(B/((zeta+B)^3)*y(1)*y(2))+(3*B/((zeta+B)^2)*y(2)*y(2))+(3*B/((zeta+B)^2)*y(2)*y(3))+2*FR(i)*(2*y(2)*y(3)+((1/(zeta+B))*y(2)*y(2))));
dydx(5) = y(6); % y(5)=theta
dydx(6) = -(1/(zeta+B))*y(6)-((Pr*shi3*B)/(Kb*(zeta+B)))*(y(1)*y(6)-A*y(2)*y(5));% k=khnf/kbf
%--------------------------------------------------------------------------
end
function res = bcs1(ya,yb)
res = [ya(1); ya(2)-1;ya(5)-1; yb(2); yb(3); yb(5)];
end
end

采纳的回答

Torsten
Torsten 2023-9-13
编辑:Torsten 2023-9-13
S_a()
function S_a
%%%%note%%%%%%%%%
%%%%%%%%%%%%%%%%% decreasing behaviour with Fr=0.4,0.8,1.2,1.6,2 %%%%%%%%%%%
FR = 0.4:0.4:2;
F2 = 0:0.01:0.15;
Pr = 6.2 ;
B = 0.9;
f1 = 0.05 ;
lambda = 0.4;
zeta=0.5;
A= 0.7;
n=3.7;
Re=0.3;
C1=765;
P1=3970;
K1=40;
%%%%%%%%%%%%%%%%%%%%%%%
C2=385; % specific heat
P2=8933; % density
K2=400; % thermal conductivity
%%%%%%%%%%% %%%%%%%%%%%%
C3=4180; % specific heat
P3=997.1; % density
K3=0.6071; % thermal conductivity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%
Z=zeta+B;
H1=P1*C1; % pho*cp nanoparticle 1
H2=P2*C2; % pho*cp nanoparticle 2
H3=P3*C3; % pho*cp base fluid
hold on
for j=1:numel(FR)
Sk = zeros(size(F2));
Fr = FR(j); %variation paramter
for i=1:numel(F2)
f2 = F2(i);
Kb = (K2 + (n - 1)*K3 - f2*(n - 1)*(K3 - K2))/(K2 + (n - 1)*K3 + f2*(K3 - K2)); % khnf/kbf
Kf=((K2 + (n - 1)*K3 - f2*(n - 1)*(K3 - K2))/(K2 + (n - 1)*K3 + f2*(K3 - K2)))*((K1 + (n - 1)*K3 - f1*(n - 1)*(K3 - K1))/(K1 + (n - 1)*K3 + f1*(K3 - K1)));
% khnf/kf
shi2=(((1 - f1)^2.5)*((1 - f2)^2.5))*(((1 - f1)*(1 - f2))+f1*(P1/P3)) + f2*(P2/P3) ;
% shi2 in velocity equation
shi3=(1 - f2)*((1 - f1) + f1*((H1)/(H3))) + f2*((H2)/(H3));
u=1/((1-f1)^2.5);
options = bvpset('RelTol',1e-6,'Stats','off') ;
solinit = bvpinit(linspace(0,5,10),[1; 0 ;0 ;0 ;0; 0]);
sol = bvp4c(@odes,@bcs,solinit,options);
Sk(i)=u*1/((Re)^0.5)*(sol.y(3,1)-sol.y(2,1)/B);
end
plot(F2,Sk,'-y')
end
hold off
%--------------------------------------------------------------------------
function dydx = odes(~,y)
% global zeta Pr psi2 psi3 k B Fr A
dydx = zeros(6,1);
dydx(1) = y(2); % y(1)=f
dydx(2) = y(3);
dydx(3) = y(4);
dydx(4) = -(2/(zeta+B)*y(4))+(1/((zeta+B)^2)*y(3))-(1/((zeta+B)^3)*y(2))+2*lambda*(y(3)+(1/(zeta+B)*y(2)))+shi2*(-((B/(zeta+B)^2)*y(1)*y(3))-((B/(zeta+B)*y(1)*y(4)))+(B/((zeta+B)^3)*y(1)*y(2))+(3*B/((zeta+B)^2)*y(2)*y(2))+(3*B/((zeta+B)^2)*y(2)*y(3))+2*Fr*(2*y(2)*y(3)+((1/(zeta+B))*y(2)*y(2))));
dydx(5) = y(6); % y(5)=theta
dydx(6) = -(1/(zeta+B))*y(6)-((Pr*shi3*B)/(Kb*(zeta+B)))*(y(1)*y(6)-A*y(2)*y(5));% k=khnf/kbf
%--------------------------------------------------------------------------
end
function res = bcs(ya,yb)
res = [ya(1); ya(2)-1;ya(5)-1; yb(2); yb(3); yb(5)];
end
end

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Engines & Motors 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by