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

2 次查看(过去 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 中查找有关 Heat and Mass Transfer 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by