Magnetized Hybrid Nanofluid Flow , MATLAB Code has some problem. Please help to Rectify. Highly Appreciated

function MHN

% Initialization of paramters

beta=1.5;

lambda=1.5;

wt=2.5;

wb=1.5;

ks1=0.5;

ks2=0.1;

a=0.5;

epsilon=0.1;

delta1=0.1;

rhos1=0.2;

rhos2=0.3;

omegas1=0.2;

omegas2=0.1;

phi1=0.5;

phi2=0.5;

rhocps1=0.3;

rhocps2=0.1;

chi=0.1;

omegaf=0.05;

rhof=997.1;

kf=0.613;

rhocpf=4179;

% Constants involve in Equation # 14

CC1=((1-phi1)^(2.5)).*((1-phi2)^(2.5)); % Hybrid to Nanofluid Constant

E1=(1/CC1); % Mau_Hnf/Mau_f (Equation # 14)

CC2=(1-phi2).*((1-phi1).*rhof + rhos1.*phi1) + rhos2.*phi2;% Hybrid to Nanofluid Constant

E2=CC2.*(1/rhof); % rho_Hnf/rho_f (Equation # 14)

DD1=omegas2.*(1+2.*phi2)+2.*omegaf.*(1-phi2); % Hybrid Constant

DD2=omegas2.*(1-phi2)+omegaf.*(2+phi2); % Hybrid

CC3=DD1/DD2; % Hybrid Constant

DD3=omegas1.*(1+2.*phi1)+2.*omegaf.*(1-phi1); % Nanofluid Constant

DD4=omegas1.*(1-phi1)+omegaf.*(2+phi1); % Nanofluid Constant

CC4=DD3/DD4; % Nanofluid Constant

E3=CC3.*CC4; % omega_Hnf/omega_f (Equation # 14)

EE1=2.*kf+ks1-2.*(kf-ks1).*phi1; % Nanofluid Constant

EE2=2.*kf+ks1+(kf-ks1).*phi1; % Nanofluid Constant

CC5=EE1/EE2; % Nanofluid Constant

EE3=2.*kf+ks2-2.*(kf-ks2).*phi2; % Hybrid Constant

EE4=2.*kf+ks2+(kf-ks2).*phi2; % Hybrid Constant

CC6=EE3/EE4; % Hybrid Constant

E4=CC5.*CC6; % k_Hnf/k_f (Equation # 14)

FF1=(1-phi2)*((1-phi1)*rhocpf+phi1*rhocps1)+phi2*rhocps2; % Hybrid to Nanofluid Constant

FF2=1/rhocpf; % Hybrid to Nanofluid Constant

E5=FF1.*FF2; % rhocp_Hnf/rhocp_f (Equation # 14)

% Initial Condition Input

sol = bvpinit(linspace(0,5,10), [1 0 0 0 0 0 0]);

% solution in structure form

sol1 = bvp4c(@bvpexam2, @bcexam2, sol);

x1 = sol1.x;

y1 = sol1.y;

plot(x1, y1(2,:));

figure (1)

hold on

value = deval(sol1,0);

vpa(value,9);

function res=bcexam2(y0, yinf)

res=[y0(1);y0(2)-1;y0(4)-1-delta1*y0(5);y0(6)-1; yinf(2);yinf(4);yinf(6)]

end

function dydx = bvpexam2(t,y)

yy1=(E3/E1)*beta*y(2)-(E2/E1)*y(1)*y(3)+(E2/E1)*y(2)^(2)

yy2 = -(chi/E4(1+epsilon*y(4)))*(E5*y(1)*y(5)+wt*y(5)*y(7)+wb*y(5)^(2)+lambda(E1*y(3)^(2)+E3*beta*y(2)^(2))+epsilon*E4*y(5)^(2))

yy3 = -a*(y(1)*y(7))-(wt/wb)*yy2

dydx= [y(2);y(3);yy1;y(5);yy2;y(7);yy3]

end

Error in MATLAB

Array indices must be positive integers or logical values.

Error in MHN/bvpexam2 (line 62)

yy2 =

-(chi/E4(1+epsilon*y(4)))*(E5*y(1)*y(5)+wt*y(5)*y(7)+wb*y(5)*y(5)+lambda(E1*y(3)*y(3)+E3*beta*y(2)*y(2))+epsilon*E4*y(5)*y(5))

Error in bvparguments (line 105)

testODE = ode(x1,y1,odeExtras{:});

Error in bvp4c (line 128)

bvparguments(solver_name,ode,bc,solinit,options,varargin);

Error in MHN (line 48)

sol1 = bvp4c(@bvpexam2, @bcexam2, sol);

[Shahid Hasnain]

Attached file has equation # 14 for understanding of parameters and system of ODEs.

[Fareeha]

@shahid is this code working well for temperature profile? As i modified it for another problem and it is not showing proper results. Thank you in advance

[Shahid Hasnain]

@Fareeha

I am uncertain about the specific issue you are referring to. Furthermore, these modifications involve not only the system of equations but also the boundary conditions and associated parameters. I trust that you will gain clarity regarding the specific inquiry you wish to make.

[Fareeha]

@Shahid I am sorry for not sharing details here. Here is the link to my query. Any help resolving the issue will be appreciated. https://www.mathworks.com/matlabcentral/answers/2007417-yy3-is-generating-a-straight-line-while-it-should-generate-a-curve-i-don-t-understand-what-is-wrong?s_tid=srchtitle

Answer 1

[Torsten]

MHN

function MHN
% Initialization of paramters
beta=1.5;
lambda=1.5;
wt=2.5;
wb=1.5;
ks1=0.5;
ks2=0.1;
a=0.5;
epsilon=0.1;
delta1=0.1;
rhos1=0.2;
rhos2=0.3;
omegas1=0.2;
omegas2=0.1;
phi1=0.5;
phi2=0.5;
rhocps1=0.3;
rhocps2=0.1;
chi=0.1;
omegaf=0.05;
rhof=997.1;
kf=0.613;
rhocpf=4179;
% Constants involve in Equation # 14
CC1=((1-phi1)^(2.5)).*((1-phi2)^(2.5));                    % Hybrid to Nanofluid Constant
E1=(1/CC1);                                                % Mau_Hnf/Mau_f (Equation # 14)
CC2=(1-phi2).*((1-phi1).*rhof + rhos1.*phi1) + rhos2.*phi2;% Hybrid to Nanofluid Constant
E2=CC2.*(1/rhof);                                          % rho_Hnf/rho_f (Equation # 14)
DD1=omegas2.*(1+2.*phi2)+2.*omegaf.*(1-phi2);              % Hybrid Constant
DD2=omegas2.*(1-phi2)+omegaf.*(2+phi2);                    % Hybrid
CC3=DD1/DD2;                                               % Hybrid Constant
DD3=omegas1.*(1+2.*phi1)+2.*omegaf.*(1-phi1);              % Nanofluid Constant
DD4=omegas1.*(1-phi1)+omegaf.*(2+phi1);                    % Nanofluid Constant
CC4=DD3/DD4;                                               % Nanofluid Constant
E3=CC3.*CC4;                                               % omega_Hnf/omega_f (Equation # 14)
EE1=2.*kf+ks1-2.*(kf-ks1).*phi1;                           % Nanofluid Constant
EE2=2.*kf+ks1+(kf-ks1).*phi1;                              % Nanofluid Constant
CC5=EE1/EE2;                                               % Nanofluid Constant
EE3=2.*kf+ks2-2.*(kf-ks2).*phi2;                           % Hybrid Constant
EE4=2.*kf+ks2+(kf-ks2).*phi2;                              % Hybrid Constant
CC6=EE3/EE4;                                               % Hybrid Constant
E4=CC5.*CC6;                                               % k_Hnf/k_f (Equation # 14)
FF1=(1-phi2)*((1-phi1)*rhocpf+phi1*rhocps1)+phi2*rhocps2;  % Hybrid to Nanofluid Constant
FF2=1/rhocpf;                                              % Hybrid to Nanofluid Constant
E5=FF1.*FF2;                                               % rhocp_Hnf/rhocp_f (Equation # 14)

% Initial Condition Input
sol = bvpinit(linspace(0,5,10), [1 0 0 0 0 0 0]);
%  solution in structure form
sol1 = bvp4c(@bvpexam2, @bcexam2, sol);
x1 = sol1.x;
y1 = sol1.y;
plot(x1, y1(2,:));
figure (1)
hold on
value = deval(sol1,0);
vpa(value,9);

    function res=bcexam2(y0, yinf)
        res=[y0(1);y0(2)-1;y0(4)-1-delta1*y0(5);y0(6)-1; yinf(2);yinf(4);yinf(6)];  
    end
    function dydx = bvpexam2(t,y)
        yy1=(E3/E1)*beta*y(2)-(E2/E1)*y(1)*y(3)+(E2/E1)*y(2)^(2);
        yy2 = -(chi/E4*(1+epsilon*y(4)))*(E5*y(1)*y(5)+wt*y(5)*y(7)+wb*y(5)^(2)+lambda*(E1*y(3)^(2)+E3*beta*y(2)^(2))+epsilon*E4*y(5)^(2));
        yy3 = -a*(y(1)*y(7))-(wt/wb)*yy2;
        dydx= [y(2);y(3);yy1;y(5);yy2;y(7);yy3] ;       
    end
end

[Shahid Hasnain]

It is working perfactly. Many many Thanks @Torsten

[Fareeha]

@Torsten can we automate this code for ploting the figure for different values of beta?

[Shahid Hasnain]

@Fareeha

I guess use for loop, like

TT=[1.1, 1.2, 1.5, 1.9, 2.3];

for i=1:length(TT)

beta= TT(i);

whole code must be here

end

[Fareeha]

@Shahid do you mean

function MHM
% Initialization of paramters
beta=1.5;
lambda=1.5;
wt=2.5;
wb=1.5;
ks1=0.5;
ks2=0.1;
a=0.5;
epsilon=0.1;
delta1=0.1;
rhos1=0.2;
rhos2=0.3;
omegas1=0.2;
omegas2=0.1;
phi1=0.5;
phi2=0.5;
rhocps1=0.3;
rhocps2=0.1;
chi=0.1;
omegaf=0.05;
rhof=997.1;
kf=0.613;
rhocpf=4179;
% Constants involve in Equation # 14
CC1=((1-phi1)^(2.5)).*((1-phi2)^(2.5)); % Hybrid to Nanofluid Constant
E1=(1/CC1); % Mau_Hnf/Mau_f (Equation # 14)
CC2=(1-phi2).*((1-phi1).*rhof + rhos1.*phi1) + rhos2.*phi2;% Hybrid to Nanofluid Constant
E2=CC2.*(1/rhof); % rho_Hnf/rho_f (Equation # 14)
DD1=omegas2.*(1+2.*phi2)+2.*omegaf.*(1-phi2); % Hybrid Constant
DD2=omegas2.*(1-phi2)+omegaf.*(2+phi2); % Hybrid
CC3=DD1/DD2; % Hybrid Constant
DD3=omegas1.*(1+2.*phi1)+2.*omegaf.*(1-phi1); % Nanofluid Constant
DD4=omegas1.*(1-phi1)+omegaf.*(2+phi1); % Nanofluid Constant
CC4=DD3/DD4; % Nanofluid Constant
E3=CC3.*CC4; % omega_Hnf/omega_f (Equation # 14)
EE1=2.*kf+ks1-2.*(kf-ks1).*phi1; % Nanofluid Constant
EE2=2.*kf+ks1+(kf-ks1).*phi1; % Nanofluid Constant
CC5=EE1/EE2; % Nanofluid Constant
EE3=2.*kf+ks2-2.*(kf-ks2).*phi2; % Hybrid Constant
EE4=2.*kf+ks2+(kf-ks2).*phi2; % Hybrid Constant
CC6=EE3/EE4; % Hybrid Constant
E4=CC5.*CC6; % k_Hnf/k_f (Equation # 14)
FF1=(1-phi2)*((1-phi1)*rhocpf+phi1*rhocps1)+phi2*rhocps2; % Hybrid to Nanofluid Constant
FF2=1/rhocpf; % Hybrid to Nanofluid Constant
E5=FF1.*FF2; % rhocp_Hnf/rhocp_f (Equation # 14)
TT=[1.1, 1.2, 1.5, 1.9, 2.3];
for i=1:length(TT)
    beta= TT(i);
    % Initial Condition Input
    sol = bvpinit(linspace(0,5,10), [1 0 0 0 0 0 0]);
    % solution in structure form
    sol1 = bvp4c(@bvpexam2, @bcexam2, sol);
    x1 = sol1.x;
    y1 = sol1.y;
    plot(x1, y1(2,:));
    figure (1)
    hold on
    value = deval(sol1,0);
    vpa(value,9);
    function res=bcexam2(y0, yinf)
        res=[y0(1);y0(2)-1;y0(4)-1-delta1*y0(5);y0(6)-1; yinf(2);yinf(4);yinf(6)]; 
    end
    function dydx = bvpexam2(t,y)
        yy1=(E3/E1)*beta*y(2)-(E2/E1)*y(1)*y(3)+(E2/E1)*y(2)^(2);
        yy2 = -(chi/E4*(1+epsilon*y(4)))*(E5*y(1)*y(5)+wt*y(5)*y(7)+wb*y(5)^(2)+lambda*(E1*y(3)^(2)+E3*beta*y(2)^(2))+epsilon*E4*y(5)^(2));
        yy3 = -a*(y(1)*y(7))-(wt/wb)*yy2;
        dydx= [y(2);y(3);yy1;y(5);yy2;y(7);yy3] ; 
    end
end

[Fareeha]

@Shahid I have tried this within function dydx block function dydx = bvpexam2(t,y) TT=[1.1, 1.2, 1.5, 1.9, 2.3]; for i=1:length(TT) beta= TT(i); yy1=(E3/E1)*beta*y(2)-(E2/E1)*y(1)*y(3)+(E2/E1)*y(2)^(2); yy2 = -(chi/E4*(1+epsilon*y(4)))*(E5*y(1)*y(5)+wt*y(5)*y(7)+wb*y(5)^(2)+lambda*(E1*y(3)^(2)+E3*beta*y(2)^(2))+epsilon*E4*y(5)^(2)); yy3 = -a*(y(1)*y(7))-(wt/wb)*yy2; dydx= [y(2);y(3);yy1;y(5);yy2;y(7);yy3] ; end end And getting this output

</matlabcentral/answers/uploaded_files/1440283/Figure%202023-07-23%2010_41_42.png> Where I am mistaken?

[Shahid Hasnain]

@Fareeha

function MHN

% Initialization of paramters

Megnatic_C=[0.5 0.9 1.3 1.5 1.9]

for i=1:length(Megnatic_C)

beta=Megnatic_C(i);

lambda=1.5;

wt=2.5;

wb=1.5;

ks1=0.5;

ks2=0.1;

a=0.5;

epsilon=0.1;

delta1=0.1;

rhos1=0.2;

rhos2=0.3;

omegas1=0.2;

omegas2=0.1;

phi1=0.5;

phi2=0.5;

rhocps1=0.3;

rhocps2=0.1;

chi=0.1;

omegaf=0.05;

rhof=997.1;

kf=0.613;

rhocpf=4179;

% Constants involve in Equation # 14

CC1=((1-phi1)^(2.5)).*((1-phi2)^(2.5)); % Hybrid to Nanofluid Constant

E1=(1/CC1); % Mau_Hnf/Mau_f (Equation # 14)

CC2=(1-phi2).*((1-phi1).*rhof + rhos1.*phi1) + rhos2.*phi2;% Hybrid to Nanofluid Constant

E2=CC2.*(1/rhof); % rho_Hnf/rho_f (Equation # 14)

DD1=omegas2.*(1+2.*phi2)+2.*omegaf.*(1-phi2); % Hybrid Constant

DD2=omegas2.*(1-phi2)+omegaf.*(2+phi2); % Hybrid

CC3=DD1/DD2; % Hybrid Constant

DD3=omegas1.*(1+2.*phi1)+2.*omegaf.*(1-phi1); % Nanofluid Constant

DD4=omegas1.*(1-phi1)+omegaf.*(2+phi1); % Nanofluid Constant

CC4=DD3/DD4; % Nanofluid Constant

E3=CC3.*CC4; % omega_Hnf/omega_f (Equation # 14)

EE1=2.*kf+ks1-2.*(kf-ks1).*phi1; % Nanofluid Constant

EE2=2.*kf+ks1+(kf-ks1).*phi1; % Nanofluid Constant

CC5=EE1/EE2; % Nanofluid Constant

EE3=2.*kf+ks2-2.*(kf-ks2).*phi2; % Hybrid Constant

EE4=2.*kf+ks2+(kf-ks2).*phi2; % Hybrid Constant

CC6=EE3/EE4; % Hybrid Constant

E4=CC5.*CC6; % k_Hnf/k_f (Equation # 14)

FF1=(1-phi2)*((1-phi1)*rhocpf+phi1*rhocps1)+phi2*rhocps2; % Hybrid to Nanofluid Constant

FF2=1/rhocpf; % Hybrid to Nanofluid Constant

E5=FF1.*FF2; % rhocp_Hnf/rhocp_f (Equation # 14)

% Initial Condition Input

sol = bvpinit(linspace(0,5,10), [1 0 0 0 0 0 0]);

% solution in structure form

sol1 = bvp4c(@bvpexam2, @bcexam2, sol);

x1 = sol1.x;

y1 = sol1.y;

plot(x1, y1(2,:));

figure (1)

hold on

value = deval(sol1,0);

vpa(value,9);

end

function res=bcexam2(y0, yinf)

res=[y0(1);y0(2)-1;y0(4)-1-delta1*y0(5);y0(6)-1; yinf(2);yinf(4);yinf(6)];

end

function dydx = bvpexam2(t,y)

yy1=(E3/E1)*beta*y(2)-(E2/E1)*y(1)*y(3)+(E2/E1)*y(2)^(2)

yy2 = -(chi/E4*(1+epsilon*y(4)))*(E5*y(1)*y(5)+wt*y(5)*y(7)+wb*y(5)^(2)+lambda*(E1*y(3)^(2)+E3*beta*y(2)^(2))+epsilon*E4*y(5)^(2))

yy3 = -a*(y(1)*y(7))-(wt/wb)*yy2

dydx= [y(2);y(3);yy1;y(5);yy2;y(7);yy3]

end

end

[Fareeha]

@Shahid Hasnain Thank you so much