sol = bvp4c (OdeBVP, OdeBC, solinit, options);

118 次查看(过去 30 天)
ne the boundary conditions
function res = OdeBc (ya, yb, A, s, B, lambda)
global A s B lambda
res= [ya(1)-s;
ya(2)-lambda-A*ya(3);
ya(4)-1-B*ya(5);
yb(2);
yb(4)];
end
% setting the initial guess for first solution
function v = OdeInit1(x,A,s,lambda)
global A s lambda
v=[s+0.56
0
0
0
0];
end
% setting the initial guess for second solution
function v1 =OdeInit2(x, A, s)
global A s
v1 = [exp(-x)
exp(-x)
-exp(-x)
-exp(-x)
-exp(-x)];
end
end
  16 个评论
Waseef
Waseef 2024-6-7
sorry sir for inconveniance, how i generate stream lines for this problem like in the image thank you.
Torsten
Torsten 2024-6-7
In your code, η is between 0 and 15 instead of 0 and 1 and ξ is not defined.

请先登录,再进行评论。

采纳的回答

Torsten
Torsten 2022-5-31
function slipflow
format long g
%Define all parameters
% Boundary layer thickness & stepsize
etaMin = 0;
etaMax1 = 15;
etaMax2 = 15; %15, 10
stepsize1 = etaMax1;
stepsize2 = etaMax2;
% Input for the parameters
A=0.6; %velocity slip
B=0.2; %thermal slip
beta=0.02; %heat gen/abs
S=2.4; %suction(2.3,2.4,2.5)
Pr=6.2; %prandtl number
lambda=-1; %stretching shrinking
a=0.01; %phil-1st nanoparticle concentration
b=0.01; %(0.01,0.05)phi2-2nd nanoparticle concentration
c=a+b; %phi-hnf concentration of hybrid nanoparticle
%%%%%%%%%%% 1st nanoparticle properties (Al2O3)%%%%%%%%%%%%
C1=765;
P1=3970;
K1=40;
B1=0.85/((10)^5);
s1=35*(10)^6; %MHD
%%%%%%%%%%% 2nd nanoparticle properties (Cu)%%%%%%%%%%%%
C2=385; %specific heat
P2=8933; %density
K2=400; %thermal conductivity
B2=1.67/((10)^5); %thermal expansion
s2=(59.6)*(10)^6; %MHD
%%%%%%%%%%% Base fluid properties %%%%%%%%%%%%
C3=4179; %specific heat
P3=997.1; %density
K3=0.613; %thermal conductivity
B3=21/((10)^5); %thermal expansion
s3=0.05; %MHD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%multiplier%%%%%%%%%%%%%%%%%%%
H1=P1*C1; %pho*cp nanoparticle 1
H2=P2*C2; %pho*cp nanoparticle 2
H3=P3*C3; %pho*cp base fluid
H4=a*H1+b*H2+(1-c)*H3; %pho*cp hybrid nanofluid
H5=a*P1+b*P2+(1-c)*P3; %pho hybrid nanofluid
H6=1/((1-c)^2.5); % mu hybrid nanofluid / mu base fluid
H7=a*(P1*B1)+b*(P2*B2)+(1-c)*(P3*B3); % thermal expansion of hybrid nanofluid
%Kn=K3*(K1+2*K3-2*a*(K3-K1))/(K1+2*K3+a*(K3-K1)); %thermal conductivity of nanofluid
Kh=(((a*K1+b*K2)/c)+2*K3+2*(a*K1+b*K2)-2*c*K3)/(((a*K1+b*K2)/c)+2*K3-(a*K1+b*K2)-2*c*K3); %khnf/kf
H8=(((a*s1+b*s2)/c)+2*s3+2*(a*s1+b*s2)-2*c*s3)/(((a*s1+b*s2)/c)+2*s3-(a*s1+b*s2)-2*c*s3); % \sigma hnf/ \sigma f
D1=(H5/P3)/H6;
D3=(H7/(P3*B3))/(H5/P3); % multiplier of boundary parameter
D2= Pr*((H4/H3)/Kh);
D4=H8/(H5/P3); %multiplier MHD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First solution %%%%%%%%%%%%%%%%%%%
options = bvpset('stats','off','RelTol',1e-10);
solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1),@(x)OdeInit1(x,A,S,lambda));
sol = bvp4c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBC(ya, yb, A, S, B, lambda), solinit, options);
eta = linspace (etaMin, etaMax1, stepsize1);
y= deval (sol,eta);
figure(1) %velocity profile
plot(sol.x,sol.y(2,:),'-')
xlabel('\eta')
ylabel('f`(\eta)')
hold on
figure(2) %temperature profile
plot(sol.x,sol.y(4,:),'-')
xlabel('\eta')
ylabel('\theta(\eta)')
hold on
% saving the out put in text file for first solution
descris =[sol.x; sol.y];
save 'sliphybrid_upper.txt' descris -ascii
% Displaying the output for first solution
fprintf('\n First solution:\n');
fprintf('f"(0)=%7.9f\n',y(3)); % reduced skin friction
fprintf('-theta(0)=%7.9f\n',-y(5)); %reduced local nusselt number
fprintf('Cfx=%7.9f\n',H6*(y(3))); % skin friction
fprintf('Nux=%7.9f\n',-Kh*y(5)); % local nusselt number
fprintf('\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% second solution %%%%%%%%%%%%%%%%%%%
options = bvpset('stats','off','RelTol',1e-10);
solinit = bvpinit (linspace (etaMin, etaMax2, stepsize2),@(x)OdeInit2(x,A,S,lambda));
sol = bvp4c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBC(ya, yb, A, S, B, lambda), solinit, options);
eta= linspace (etaMin, etaMax2, stepsize2);
y = deval (sol,eta);
figure(1) %velocity profile
plot(sol.x,sol.y(2,:),'--')
xlabel('\eta')
ylabel('f`(\eta)')
hold on
figure(2) %temperature profile
plot(sol.x,sol.y(4,:),'--')
xlabel('\eta')
ylabel('\theta(\eta)')
hold on
% saving the out put in text file for second solution
descris=[sol.x; sol.y];
save 'sliphybrid_lower.txt'descris -ascii
% Displaying the output for first solution
fprintf('\nSecond solution:\n');
fprintf('f"(0)=%7.9f\n',y(3)); % reduced skin friction
fprintf('-theta(0)=%7.9f\n',-y(5)); %reduced local nusselt number
fprintf('Cfx=%7.9f\n',H6*(y(3))); % skin friction
fprintf('Nux=%7.9f\n',-Kh*y(5)); % local nusselt number
fprintf('\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
% Define the ODE function
function f = OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta)
f =[y(2);y(3);D1*(2*(y(2)*y(2))-y(1)*y(3));y(5);(Pr/Kh)*((-H4/H3)*(y(1)*y(5)-y(2)*y(4))-beta*y(4))];
end
% Define the boundary conditions
function res = OdeBc (ya, yb, A, S, B, lambda)
res= [ya(1)-S;ya(2)-lambda-A*ya(3);ya(4)-1-B*ya(5);yb(2);yb(4)];
end
% setting the initial guess for first solution
function v = OdeInit1(x,A,S,lambda)
v=[S+0.56;0;0;0;0];
end
% setting the initial guess for second solution
function v1 =OdeInit2(x, A, S,lambda)
v1 = [exp(-x);exp(-x);-exp(-x);-exp(-x);-exp(-x)];
end
  1 个评论
Torsten
Torsten 2022-5-31
编辑:Torsten 2022-5-31
Rename the function
OdeBc
in
OdeBC
For MATLAB, small and capital letters define different things.
Same with the letter "s" you used in your original code. I hope you always meant your "S" - I replaced it where it was needed.

请先登录,再进行评论。

更多回答(2 个)

Farooq Aamir
Farooq Aamir 2023-9-1
编辑:Torsten 2023-9-1
This working now.
slipflow()
First solution: f"(0)=0.954395347 -theta(0)=3.658786212 Cfx=1.003836827 Nux=3.881552117
Second solution: f"(0)=0.075669302 -theta(0)=3.617987101 Cfx=0.079589273 Nux=3.838268944
function slipflow
format long g
%Define all parameters
% Boundary layer thickness & stepsize
etaMin = 0;
etaMax1 = 15;
etaMax2 = 15; %15, 10
stepsize1 = etaMax1;
stepsize2 = etaMax2;
% Input for the parameters
A=0.6; %velocity slip
B=0.2; %thermal slip
beta=0.02; %heat gen/abs
S=2.4; %suction(2.3,2.4,2.5)
Pr=6.2; %prandtl number
lambda=-1; %stretching shrinking
a=0.01; %phil-1st nanoparticle concentration
b=0.01; %(0.01,0.05)phi2-2nd nanoparticle concentration
c=a+b; %phi-hnf concentration of hybrid nanoparticle
%%%%%%%%%%% 1st nanoparticle properties (Al2O3)%%%%%%%%%%%%
C1=765;
P1=3970;
K1=40;
B1=0.85/((10)^5);
s1=35*(10)^6; %MHD
%%%%%%%%%%% 2nd nanoparticle properties (Cu)%%%%%%%%%%%%
C2=385; %specific heat
P2=8933; %density
K2=400; %thermal conductivity
B2=1.67/((10)^5); %thermal expansion
s2=(59.6)*(10)^6; %MHD
%%%%%%%%%%% Base fluid properties %%%%%%%%%%%%
C3=4179; %specific heat
P3=997.1; %density
K3=0.613; %thermal conductivity
B3=21/((10)^5); %thermal expansion
s3=0.05; %MHD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%multiplier%%%%%%%%%%%%%%%%%%%
H1=P1*C1; %pho*cp nanoparticle 1
H2=P2*C2; %pho*cp nanoparticle 2
H3=P3*C3; %pho*cp base fluid
H4=a*H1+b*H2+(1-c)*H3; %pho*cp hybrid nanofluid
H5=a*P1+b*P2+(1-c)*P3; %pho hybrid nanofluid
H6=1/((1-c)^2.5); % mu hybrid nanofluid / mu base fluid
H7=a*(P1*B1)+b*(P2*B2)+(1-c)*(P3*B3); % thermal expansion of hybrid nanofluid
%Kn=K3*(K1+2*K3-2*a*(K3-K1))/(K1+2*K3+a*(K3-K1)); %thermal conductivity of nanofluid
Kh=(((a*K1+b*K2)/c)+2*K3+2*(a*K1+b*K2)-2*c*K3)/(((a*K1+b*K2)/c)+2*K3-(a*K1+b*K2)-2*c*K3); %khnf/kf
H8=(((a*s1+b*s2)/c)+2*s3+2*(a*s1+b*s2)-2*c*s3)/(((a*s1+b*s2)/c)+2*s3-(a*s1+b*s2)-2*c*s3); % \sigma hnf/ \sigma f
D1=(H5/P3)/H6;
D3=(H7/(P3*B3))/(H5/P3); % multiplier of boundary parameter
D2= Pr*((H4/H3)/Kh);
D4=H8/(H5/P3); %multiplier MHD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First solution %%%%%%%%%%%%%%%%%%%
options = bvpset('stats','off','RelTol',1e-10);
solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1),@(x)OdeInit1(x,A,S,lambda));
sol = bvp4c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBc(ya, yb, A, S, B, lambda), solinit, options);
eta = linspace (etaMin, etaMax1, stepsize1);
y= deval (sol,eta);
figure(1) %velocity profile
plot(sol.x,sol.y(2,:),'-')
xlabel('\eta')
ylabel('f`(\eta)')
hold on
figure(2) %temperature profile
plot(sol.x,sol.y(4,:),'-')
xlabel('\eta')
ylabel('\theta(\eta)')
hold on
% saving the out put in text file for first solution
descris =[sol.x; sol.y];
%save 'sliphybrid_upper.txt' descris -ascii
% Displaying the output for first solution
fprintf('\n First solution:\n');
fprintf('f"(0)=%7.9f\n',y(3)); % reduced skin friction
fprintf('-theta(0)=%7.9f\n',-y(5)); %reduced local nusselt number
fprintf('Cfx=%7.9f\n',H6*(y(3))); % skin friction
fprintf('Nux=%7.9f\n',-Kh*y(5)); % local nusselt number
fprintf('\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% second solution %%%%%%%%%%%%%%%%%%%
options = bvpset('stats','off','RelTol',1e-10);
solinit = bvpinit (linspace (etaMin, etaMax2, stepsize2),@(x)OdeInit2(x,A,S,lambda));
sol = bvp4c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBc(ya, yb, A, S, B, lambda), solinit, options);
eta= linspace (etaMin, etaMax2, stepsize2);
y = deval (sol,eta);
figure(1) %velocity profile
plot(sol.x,sol.y(2,:),'--')
xlabel('\eta')
ylabel('f`(\eta)')
hold on
figure(2) %temperature profile
plot(sol.x,sol.y(4,:),'--')
xlabel('\eta')
ylabel('\theta(\eta)')
hold on
% saving the out put in text file for second solution
descris=[sol.x; sol.y];
%save 'sliphybrid_lower.txt'descris -ascii
% Displaying the output for first solution
fprintf('\nSecond solution:\n');
fprintf('f"(0)=%7.9f\n',y(3)); % reduced skin friction
fprintf('-theta(0)=%7.9f\n',-y(5)); %reduced local nusselt number
fprintf('Cfx=%7.9f\n',H6*(y(3))); % skin friction
fprintf('Nux=%7.9f\n',-Kh*y(5)); % local nusselt number
fprintf('\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
% Define the ODE function
function f = OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta)
f =[y(2);y(3);D1*(2*(y(2)*y(2))-y(1)*y(3));y(5);(Pr/Kh)*((-H4/H3)*(y(1)*y(5)-y(2)*y(4))-beta*y(4))];
end
% Define the boundary conditions
function res = OdeBc(ya, yb, A, S, B, lambda)
res= [ya(1)-S;ya(2)-lambda-A*ya(3);ya(4)-1-B*ya(5);yb(2);yb(4)];
end
% setting the initial guess for first solution
function v = OdeInit1(x,A,S,lambda)
v=[S+0.56;0;0;0;0];
end
% setting the initial guess for second solution
function v1 =OdeInit2(x, A, S,lambda)
v1 = [exp(-x);exp(-x);-exp(-x);-exp(-x);-exp(-x)];
end
  9 个评论
Ramanuja
Ramanuja 2024-3-25
Thanks Farooq Aamir sir,
Thanks Farooq Aamir sir,
Thanks Farooq Aamir sir,
Yasir
Yasir 2024-6-27,11:56
Hello sir, What changes should i make in the code to plot the graphs of skin friction,Nusselt number and sherword number.
how can i get the different [oints to plot them

请先登录,再进行评论。


Waseef
Waseef 2024-6-30,13:39
编辑:Walter Roberson 2024-7-4,4:58
sir how i define entropy in this code
this the entropy "NN=y(6)*y(6))+C1*y(7)+D*(y(5)*y(5)+y(2)*y(2))+E*(y(1)*y(8)+y(4)*y(8))"
and this is the code
Skinforbydirectional()
f"(0)=-1.145548435 g"(0)=-0.114554844 -theta(0)=1.029922860
function Skinforbydirectional
format long g
% Boundary layer thickness & stepsize
global A Pr aa pm Phi R M pm Q sigmaf sigmas sigmanf n ks kf Rhos Rhof
global Cps Cpf tt kk Ec spn phi1 phi5 phi2 Lambda A B C st A2 A3 A4 A5 A6 total
etaMin = 0;
etaMax = 10;
stepsize = etaMax;%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Define all parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Phi=0.04; %input ('Input the value of phi = ');
R =0.2; %input ('Input the value of Radiation = ');
M= 0.01; %input('magentic parameter M =');
pm=0.2; %input('porosity =');
Q=0.2; %input('Heat sourse parameter');
%alpha=0.3;
aa=0.5;
%n=3; %input ('Input the value of n = ');
Ec =0; %input ('Input Eckret for velociy exponent parameter = ');
Pr = 6.8; % input ('Input the Prandtl number = ');
st=0.1;
%aa=10;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigmaf=0.05; sigmas=1000000;
Ks= 400; Kf= 0.613;
Rhos= 8933; Rhof= 997.1;Cps= 385; Cpf= 4179; tt=Rhos*Cps;
kk=Rhof*Cpf;
spn=(3*Phi*((sigmas/sigmaf)-1))/(((sigmas/sigmaf)+2)-((sigmas/sigmaf)-1)*Phi);
sigmanf=sigmaf*(1+spn);
% Lambda = 1:1:10;
%for ii = 1:numel(Lambda) %stretching shrinking
% aa = Lambda(ii);
total=(sigmas/sigmaf);
phi3=1+(3*(total-1)*Phi/((total+2)-(total-1)*Phi));
%phi4=(ks+(n-1)*kf+(n-1)*(ks-kf)*Phi)/(ks+(n-1)*kf+(kf-ks)*Phi);
phi4=((1-Phi)+2*Phi*Ks/(Ks-Kf)*log((Ks+Kf)/2*Kf))/((1-Phi)+2*Phi*Kf/(Ks-Kf)*log((Ks+Kf)/2*Kf));
%phi5=1-Phi+Phi*(Cps/Cpf);
phi1=(1-Phi)^2.5;
phi2=1-Phi+Phi*(Rhos/Rhof);
phi5=1-Phi+Phi*(tt/kk);
A = ((phi1 * M * (sigmanf / sigmaf)) + pm) * (2 / (aa + 1));
B = phi1 * phi2 * (2 * aa / (aa + 1));
C = phi1 * phi2;
A1 = phi4 + R;
B1 = Ec * Pr / phi1;
C1 = Pr * Q * (2 / (aa + 1));
E = Pr * phi5;
D = Pr * Ec * M * (sigmanf / sigmaf) * (2 / (aa + 1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First solution %%%%%%%%%%%%%%%%%%%
options = bvpset('stats','off','RelTol',1e-10);
solinit = bvpinit (linspace (etaMin, etaMax, stepsize),@(x)OdeInit1);
sol=bvp4c(@OdeBVP, @OdeBC, solinit);
%sol = bvp5c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBc(ya, yb, A, S, B, lambda), solinit, options);
eta = linspace (etaMin, etaMax, stepsize);
y= deval (sol,eta);
% To_Plot1(ii) = (1/phi1)*(sqrt(2*(aa+1)))*sol.y(5,1);
% To_Plot2(ii) = -(phi4+R)*(sqrt((aa+1)/2))*sol.y(8,1);
% fprintf('y(3) at etaMax = %f\n', y(3, end));
% fprintf('y(8) at etaMax = %f\n', y(8, end));
%end
figure(1) %velocity profile
plot(sol.x,sol.y(2,:),'-')
xlabel('\eta')
ylabel('f`(\eta)')
hold on
figure(2) %velocity profile
plot(sol.x,sol.y(5,:),'-')
xlabel('\eta')
ylabel('q`(\eta)')
hold on
% figure(3) %velocity profile
% plot(Lambda,To_Plot1,'LineWidth',2)
% xlabel('a')
% ylabel('f^\prime^\prime(0)')
% grid on
% hold on
% figure(4) %temperature profile
% hold on
% grid on
% plot(Lambda,To_Plot2,'LineWidth',2)
% xlabel('a')
% ylabel('\theta^\prime(0)')
% %xlim([0 2])
% figure(1) %velocity profile
% plot(Lambda,To_Plot1,'-')
% xlabel('\lambda')
% ylabel('f^\prime^\prime(0)')
% xlim([1,2])
% figure(2) %temperature profile
% plot(Lambda,To_Plot2,'-')
% xlabel('\lambda')
% ylabel('\theta^\prime(0)')
% xlim([1 2])
% Define the ODE function
fprintf('f"(0)=%7.9f\n',y(3)); % reduced skin friction
fprintf('g"(0)=%7.9f\n',y(6)); % reduced skin friction
fprintf('-theta(0)=%7.9f\n',-y(8)); %reduced local nusselt number
function f = OdeBVP(~,y)
f =[ y(2); y(3);A*y(2)+B*(y(2)*y(2)+y(5)*y(2))-C*(y(1)*y(3)+y(4)*y(3));
y(5);y(6); A*y(5)+B*(y(5)*y(2)+y(5)*y(5))-C*(y(1)*y(6)+y(4)*y(6));
y(8);(-1/A1)*(B1*(y(3)*y(3)+y(6)*y(6))+C1*y(7)+D*(y(5)*y(5)+y(2)*y(2))+E*(y(1)*y(8)+y(4)*y(8)))];
end
% Define the boundary conditions
function res = OdeBC(ya, yb)
res= [ya(1);
ya(2)-1;
ya(4);
ya(5)-st;
ya(7)-1;
yb(2);
yb(5);
yb(7)];
end
% setting the initial guess for first solution
function v = OdeInit1(~)
v=[0.9
0.1
0.1
-0.1
0.1
0
-0.1
01];
end
end

类别

Help CenterFile Exchange 中查找有关 Applications 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by