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

6 Ansichten (letzte 30 Tage)
T K
T K am 10 Feb. 2023
Kommentiert: Voss am 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

Melden Sie sich an, um diese Frage zu beantworten.

Antworten (1)

Voss
Voss am 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 Kommentare
T K
T K am 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 am 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

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by