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

7 views (last 30 days)
T K
T K on 10 Feb 2023
Commented: Voss on 11 Feb 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

Sign in to answer this question.

Answers (1)

Voss
Voss on 10 Feb 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 Comments
T K
T K on 11 Feb 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 on 11 Feb 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

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!