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

조회 수: 1 (최근 30일)
farooq omer
farooq omer 2023년 9월 12일
댓글: Torsten 2023년 9월 12일
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

채택된 답변

Torsten
Torsten 2023년 9월 12일
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
Torsten
Torsten 2023년 9월 12일
편집: Torsten 2023년 9월 12일
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 2023년 9월 12일
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.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Fluid Dynamics에 대해 자세히 알아보기

제품

Community Treasure Hunt

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

Start Hunting!

Translated by