Borrar filtros
Borrar filtros

Please help me to connect between three points with plot (draw line)

6 visualizaciones (últimos 30 días)
T K
T K el 10 de Feb. de 2023
Comentada: Voss el 11 de Feb. de 2023
I want to connect the three points (draw line) shown in the picture in the plot command.
function sol= proj
clc;clf;clear;
global lamda
%Relation of base fluid
rhof=1;
kf=0.613*10^5;
cpf=4179*10^4;
muf=10^-3*10;
sigf=0.05*10^-8;
alfaf=kf/(rhof*cpf);
%FE3O4
ph1=0.01;
rho1=5180*10^-3;
cp1=670*10^4;
k1=9.7*10^5;
sig1=0.74*10^-2;
%copper
ph2=0.01;
rho2=8933*10^-3;
cp2=385*10^4;
k2=401*10^5;
sig2=5.96*10^-1;
%Relation of hyprid
m=5.7;
kh=kf*((k1+(m-1)*kf-(m-1)*ph1*(kf-k1))/((k1+(m-1)*kf+ph1*(kf-k1))))*((k2+(m-1)*kf-(m-1)*ph2*(kf-k2))/((k2+(m-1)*kf+ph2*(kf-k2))));
muh= muf/((1-ph1)^2.5*(1-ph2)^2.5);
rhoh=rhof*(1-ph2)*((1-ph1)+ph1*(rho1/rhof))+ph2*rho2;
v1 =rhof*cpf*(1-ph2)*((1-ph1)+ph1*((rho1*cp1)/(rho2*cp2)))+ph2*(rho2*cp2);
sigh=sigf+(3*((ph1*sig1+ph2*sig2)-sigf*(ph1+ph2))/(((ph1*sig1+ph2*sig2)/(sigf*(ph1+ph2)))+2-((ph1*sig1+ph2*sig2)/sigf)+(ph1+ph2)));
alfah=kh/v1;
myLegend1 = {};
rr = [0.2 0.3 0.4];
for i =1:numel(rr)
Nt = rr(i);s2=1;
qq=[13 14 15];
Prf=qq(i);
an=0.1;Lb=1;Pe=0.1;
p=-0.5;q=-0.5;M=1; L=(muf/rhof);L1=L^(p);L1=L^(p);
delta=s2/(L1);
Nb=1; Nr=1;sc=0.45;gamma=pi/3;
a=25;b=0.1;v=1;u=1;
s1=1;s2=1;
s3=sqrt(a/(muf/rhof));
Tw=273+50;Ti=273+27;deltaT=Tw-Ti;
Lf=rhof*kf;
y0 = [1,0,1,0,0,1,0,1,0,1,0];
options =bvpset('stats','on','RelTol',1e-4);
m = linspace(0,1.5);
solinit = bvpinit(m,y0);
sol= bvp4c(@projfun,@projbc,solinit,options);
figure(1)
plot(qq(i),-(sol.y(11,1)),'s')
grid on,hold on
myLegend1{i}=['Nt = ',num2str(rr(i))];
grid on,hold on
i=i+1;
end
figure(1)
legend(myLegend1)
hold on
function dy= projfun(~,y)
dy= zeros(11,1);
% alignComments
E = y(1);
dE = y(2);
F = y(3);
dF= y(4);
W = y(5);
t = y(6);
dt = y(7);
phi = y(8);
dphi = y(9);
z=y(10);
dz=y(11);
dy(1) = dE;
dy(2) = (rhoh/muh)*((((a*u)/L1^(2)))*E^2+(1/L1)*W*dE+((sigh/sigf)/(rhoh/rhof))*(1/L1^2)*M*E*sin(gamma)*sin(gamma));
dy(3) = dF;
dy(4) = (rhoh/muh)*((((b*v)/L1^(2)))*F^2+(1/L1)*W*dF+((sigh/sigf)/(rhoh/rhof))*(1/L1^2)*M*F*sin(gamma)*sin(gamma));
dy(5) = -(1/L1)*(u*a*E+b*v*F);
dy(6) = dt;
dy(7) =(Prf*(rhof/muf))*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2);
dy(8)= dphi;
dy(9)=(sc/L^(p+1))*W*dphi-(s1/s2)*(Nt/Nb)*(((Prf*(rhof/muf)))*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2));
dy(10)=dz;
dy(11)=(Lb/((muf/rhof)*L1))*W*dz+Pe*delta*dz*dphi+(Pe*delta*z+((Pe*delta*an)/(s3*(muf/rhof)^q)))*((sc/L^(p+1))*W*dphi-(s1/s2)*(Nt/Nb)*(((Prf*(rhof/muf)))*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2)));
end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
ya(3)-1;
ya(5)-0;
ya(6)-1;
ya(8)-1;ya(10)-1;
yb(1);
yb(3);
yb(6);
yb(8);
yb(10)];
end

Iniciar sesión para responder a esta pregunta.

Respuestas (1)

Voss
Voss el 10 de Feb. de 2023
Ran in:
Inside the loop, store the values necessary to plot the line after the loop. See the y_store variable below.
proj
The solution was obtained on a mesh of 282 points. The maximum residual is 7.317e-05. There were 19109 calls to the ODE function. There were 253 calls to the BC function. The solution was obtained on a mesh of 283 points. The maximum residual is 7.783e-05. There were 19152 calls to the ODE function. There were 253 calls to the BC function. The solution was obtained on a mesh of 282 points. The maximum residual is 8.303e-05. There were 19108 calls to the ODE function. There were 253 calls to the BC function.
ans = struct with fields:
solver: 'bvp4c' x: [0 0.0076 0.0152 0.0303 0.0455 0.0530 0.0606 0.0682 0.0758 0.0833 0.0909 0.0985 0.1061 0.1136 0.1212 0.1288 0.1364 0.1439 0.1515 0.1591 0.1667 0.1742 0.1818 0.1894 0.1970 0.2045 0.2121 0.2172 0.2222 0.2273 0.2323 0.2374 0.2424 0.2475 … ] y: [11×282 double] yp: [11×282 double] stats: [1×1 struct]
function sol= proj
clc;clf;clear;
global lamda
%Relation of base fluid
rhof=1;
kf=0.613*10^5;
cpf=4179*10^4;
muf=10^-3*10;
sigf=0.05*10^-8;
alfaf=kf/(rhof*cpf);
%FE3O4
ph1=0.01;
rho1=5180*10^-3;
cp1=670*10^4;
k1=9.7*10^5;
sig1=0.74*10^-2;
%copper
ph2=0.01;
rho2=8933*10^-3;
cp2=385*10^4;
k2=401*10^5;
sig2=5.96*10^-1;
%Relation of hyprid
m=5.7;
kh=kf*((k1+(m-1)*kf-(m-1)*ph1*(kf-k1))/((k1+(m-1)*kf+ph1*(kf-k1))))*((k2+(m-1)*kf-(m-1)*ph2*(kf-k2))/((k2+(m-1)*kf+ph2*(kf-k2))));
muh= muf/((1-ph1)^2.5*(1-ph2)^2.5);
rhoh=rhof*(1-ph2)*((1-ph1)+ph1*(rho1/rhof))+ph2*rho2;
v1 =rhof*cpf*(1-ph2)*((1-ph1)+ph1*((rho1*cp1)/(rho2*cp2)))+ph2*(rho2*cp2);
sigh=sigf+(3*((ph1*sig1+ph2*sig2)-sigf*(ph1+ph2))/(((ph1*sig1+ph2*sig2)/(sigf*(ph1+ph2)))+2-((ph1*sig1+ph2*sig2)/sigf)+(ph1+ph2)));
alfah=kh/v1;
myLegend1 = {};
rr = [0.2 0.3 0.4];
y_store = zeros(1,numel(rr));
for i =1:numel(rr)
Nt = rr(i);s2=1;
qq=[13 14 15];
Prf=qq(i);
an=0.1;Lb=1;Pe=0.1;
p=-0.5;q=-0.5;M=1; L=(muf/rhof);L1=L^(p);L1=L^(p);
delta=s2/(L1);
Nb=1; Nr=1;sc=0.45;gamma=pi/3;
a=25;b=0.1;v=1;u=1;
s1=1;s2=1;
s3=sqrt(a/(muf/rhof));
Tw=273+50;Ti=273+27;deltaT=Tw-Ti;
Lf=rhof*kf;
y0 = [1,0,1,0,0,1,0,1,0,1,0];
options =bvpset('stats','on','RelTol',1e-4);
m = linspace(0,1.5);
solinit = bvpinit(m,y0);
sol= bvp4c(@projfun,@projbc,solinit,options);
figure(1)
y_store(i) = -(sol.y(11,1));
plot(qq(i),-(sol.y(11,1)),'s')
grid on,hold on
myLegend1{i}=['Nt = ',num2str(rr(i))];
grid on,hold on
i=i+1;
end
figure(1)
plot(qq,y_store,'k')
legend(myLegend1)
hold on
function dy= projfun(~,y)
dy= zeros(11,1);
% alignComments
E = y(1);
dE = y(2);
F = y(3);
dF= y(4);
W = y(5);
t = y(6);
dt = y(7);
phi = y(8);
dphi = y(9);
z=y(10);
dz=y(11);
dy(1) = dE;
dy(2) = (rhoh/muh)*((((a*u)/L1^(2)))*E^2+(1/L1)*W*dE+((sigh/sigf)/(rhoh/rhof))*(1/L1^2)*M*E*sin(gamma)*sin(gamma));
dy(3) = dF;
dy(4) = (rhoh/muh)*((((b*v)/L1^(2)))*F^2+(1/L1)*W*dF+((sigh/sigf)/(rhoh/rhof))*(1/L1^2)*M*F*sin(gamma)*sin(gamma));
dy(5) = -(1/L1)*(u*a*E+b*v*F);
dy(6) = dt;
dy(7) =(Prf*(rhof/muf))*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2);
dy(8)= dphi;
dy(9)=(sc/L^(p+1))*W*dphi-(s1/s2)*(Nt/Nb)*(((Prf*(rhof/muf)))*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2));
dy(10)=dz;
dy(11)=(Lb/((muf/rhof)*L1))*W*dz+Pe*delta*dz*dphi+(Pe*delta*z+((Pe*delta*an)/(s3*(muf/rhof)^q)))*((sc/L^(p+1))*W*dphi-(s1/s2)*(Nt/Nb)*(((Prf*(rhof/muf)))*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2)));
end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
ya(3)-1;
ya(5)-0;
ya(6)-1;
ya(8)-1;ya(10)-1;
yb(1);
yb(3);
yb(6);
yb(8);
yb(10)];
end
  2 comentarios
T K
T K el 11 de Feb. de 2023
Already done, thank you very much. Finally, I would like to draw two lines at two values ​​of Nt, as shown in the image.
Voss
Voss el 11 de Feb. de 2023
You're welcome!
Where do the x- and y-coordinates of the points defining those two lines come from?
Once you have the coordinates, you can plot the lines.
https://www.mathworks.com/help/matlab/ref/plot.html

Iniciar sesión para comentar.

Categorías

Más información sobre 2-D and 3-D Plots en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by