for loop in bvp4c i get circle in plotting not solid line why?

1 visualización (últimos 30 días)
farooq omer
farooq omer el 12 de Sept. de 2023
Comentada: Torsten el 12 de Sept. de 2023
function sss
global Fr Pr f2 f1 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 = 0.7;
Pr = 6.2 ;
B = 0.9;
for f2=0:0.05:0.15
f1 = 0.05 ;
lambda = 0.4;
zeta=5;
A= 0.4;
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=u*1/((Re)^0.5)*(sol.y(3,1)-sol.y(2,1)/B);
hold on
plot(f2 ,Sk,'b-')
hold off
end
%--------------------------------------------------------------------------
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

Respuesta aceptada

Torsten
Torsten el 12 de Sept. de 2023
sss()
F2 = 1×4
0 0.0500 0.1000 0.1500
function sss
global Fr Pr f2 f1 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 = 0.7;
Pr = 6.2 ;
B = 0.9;
F2 = 0:0.05:0.15
for i=1:numel(F2)
f2 = F2(i);
f1 = 0.05 ;
lambda = 0.4;
zeta=5;
A= 0.4;
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)
%--------------------------------------------------------------------------
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
  4 comentarios
Torsten
Torsten el 12 de Sept. de 2023
Editada: Torsten el 12 de Sept. de 2023
Make sss a script instead of a function. Make odes and bcs usual functions instead of nested functions. Add the necessary global variables in the functions. Save all functions in separate files with their respective names (bcs.m, odes.m). Load the script (sss.m) into your MATLAB editor. Press the RUN button.
Torsten
Torsten el 12 de Sept. de 2023
how to draw lines for A=0.4,0.6,0.8 using a loop in this code
In the same way as you used a loop for f2.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Fluid Dynamics en Help Center y File Exchange.

Productos

Community Treasure Hunt

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

Start Hunting!

Translated by