My code is running so long and never gives solution.
3 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
clear;
close all;
g=32.2; nu=1.21e-5; %ft^2/s
QC=6;
QD=8;
QE=11;
QA=QC+QD+QE;
%Roughness for concrete e=0.001 ft
%rr=e/D
L1=600; D1=1.5; rr1=0.00067; Re1=1.7e6; Q1=11;
L2=600; D2=1.5; rr2=0.00067; Re2=1.7e6; Q2=14;
L3=800; D3=1.25; rr3=0.0008; Re3=1.7e6; Q3=7;
L4=800; D4=1.25; rr4=0.0008; Re4=1.7e6; Q4=7;
L5=400; D5=1; rr5=0.0001; Re5=1.7e6; Q5=4;
L6=400; D6=1; rr6=0.0001; Re6=1.7e6; Q6=1;
L7=800; D7=1.25; rr7=0.0008; Re7=1.7e6; Q7=5;
L8=400; D8=1; rr8=0.0001; Re8=1.7e6; Q8=1;
L9=400; D9=1; rr9=0.0001; Re9=1.7e6; Q9=4;
L10=400; D10=1; rr10=0.0001; Re10=1.7e6; Q10=4;
fF1=FrictionFactor(rr1, Re1);
fF2=FrictionFactor(rr2, Re2);
fF3=FrictionFactor(rr3, Re3);
fF4=FrictionFactor(rr4, Re4);
fF5=FrictionFactor(rr5, Re5);
fF6=FrictionFactor(rr6, Re6);
fF7=FrictionFactor(rr7, Re7);
fF8=FrictionFactor(rr8, Re8);
fF9=FrictionFactor(rr9, Re9);
fF10=FrictionFactor(rr10, Re10);
QV=[Q1; Q2; Q3; Q4; Q5; Q6; Q7; Q8; Q9; QA];
fFV=[fF1, fF2, fF3, fF4, fF5, fF6, fF7, fF8, fF9, fF10]';
eF=1;
cnf=0;
tolQ=QA*1e-6;
while eF>1e-6
cnf=cnf+1;
eQ=tolQ+1;
K1=fF1*L1*8/g/D1^5/pi^2;
K2=fF2*L2*8/g/D2^5/pi^2;
K3=fF3*L3*8/g/D3^5/pi^2;
K4=fF4*L4*8/g/D4^5/pi^2;
K5=fF5*L5*8/g/D5^5/pi^2;
K6=fF6*L6*8/g/D6^5/pi^2;
K7=fF7*L7*8/g/D7^5/pi^2;
K8=fF8*L8*8/g/D8^5/pi^2;
K9=fF9*L9*8/g/D9^5/pi^2;
K10=fF10*L10*8/g/D10^5/pi^2;
cnQ=0;
while eQ>tolQ
cnQ=cnQ+1;
F1=QV(1)+QV(2)-QA;
F2=-QV(1)+QV(3)+QV(4);
F3=QC+QV(6)-QV(2)-QV(5);
F4=QD+QV(8)-QV(4);
F5=QE-QV(6)-QV(10);
F6=QV(7)+QV(5)-QV(3);
F7=QV(9)+QV(10)-QV(8);
F8=K1*QV(1)*abs(QV(1))+K3*QV(3)*abs(QV(3))+K5*QV(5)*abs(QV(5))-K2*QV(2)*abs(QV(2));
F9=K4*QV(4)*abs(QV(4))+K8*QV(8)*abs(QV(8))+K9*QV(9)*abs(QV(9))-K7*QV(7)*abs(QV(7))-K3*QV(3)*abs(QV(3));
F10=K7*QV(7)*abs(QV(7))-K9*QV(9)*abs(QV(9))+K10*QV(10)*abs(QV(10))-K6*QV(6)*abs(QV(6))-K5*QV(5)*abs(QV(5));
Jac=[1, 1, 0, 0, 0, 0, 0, 0, 0, 0;
-1, 0, 1, 1, 0, 0, 0, 0, 0, 0;
0, -1, 0, 0, -1, 1, 0, 0, 0, 0;
0, 0, 0, -1, 0, 0, 0, 1, 0, 0;
0, 0, 0, 0, 0, -1, 0, 0, 0, -1;
0, 0, -1, 0, 1, 0, 1, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, -1, 1, 1;
2*K1*abs(QV(1)), -2*K2*abs(QV(2)), 2*K3*abs(QV(3)), 0, 2*K5*abs(QV(5)), 0, 0, 0, 0, 0;
0, 0, 0, 2*K4*abs(QV(4)), 0, 0, -2*K7*abs(QV(7)), 2*K8*abs(QV(8)), 2*K9*abs(QV(9)), 0;
0, 0, 0, 0, -2*K5*abs(QV(5)), -2*K6*abs(QV(6)), 2*K7*abs(QV(7)), 0, -2*K9*abs(QV(9)), 2*K10*abs(QV(10))];
b=[-F1; -F2; -F3; -F4; -F5; -F6; -F7; -F8; -F9; -F10];
dQV=Jac\b;
QVn=QV+dQV;
eQ=max(abs(QVn-QV));
tolQ=mean(abs(QV))/1000;
QV=QVn;
end
Re1=4*QV(1)/nu/D1/pi;
Re2=4*QV(2)/nu/D2/pi;
Re3=4*QV(3)/nu/D3/pi;
Re4=4*QV(4)/nu/D4/pi;
Re5=4*QV(5)/nu/D5/pi;
Re6=4*QV(6)/nu/D6/pi;
Re7=4*QV(7)/nu/D7/pi;
Re8=4*QV(8)/nu/D8/pi;
Re9=4*QV(9)/nu/D9/pi;
Re10=4*QV(10)/nu/D10/pi;
fF1n=FrictionFactor(rr1, Re1);
fF2n=FrictionFactor(rr2, Re2);
fF3n=FrictionFactor(rr3, Re3);
fF4n=FrictionFactor(rr4, Re4);
fF5n=FrictionFactor(rr5, Re5);
fF6n=FrictionFactor(rr6, Re6);
fF7n=FrictionFactor(rr7, Re7);
fF8n=FrictionFactor(rr8, Re8);
fF9n=FrictionFactor(rr9, Re9);
fF10n=FrictionFactor(rr10, Re10);
eF=max(abs([fF1, fF2, fF3, fF4, fF5, fF6, fF7, fF8, fF9, fF10]-[fF1n, fF2n, fF3n, fF4n, fF5n, fF6n, fF7n, fF8n, fF9n, fF10n]));
fF1=fF1n;
fF2=fF2n;
fF3=fF3n;
fF4=fF4n;
fF5=fF5n;
fF6=fF6n;
fF7=fF7n;
fF8=fF8n;
fF9=fF9n;
fF10=fF10n;
fFV=[fF1, fF2, fF3, fF4, fF5, fF6,fF7, fF8, fF9, fF10]';
end
function [f, e]=FrictionFactor(rr, Re)
e=1; f=0.01;
while e>1e-6
fn=(-0.5/log10(rr/3.7+2.51/Re/sqrt(f)))^2;
e=abs(fn-f);
f=fn;
end
end
3 Commenti
Devineni Aslesha
il 26 Ago 2019
It is mentioned in the comment that fF1=fF1n+eF; makes more sense for the issue while finding new QV values. However, this equation is not there in the provided code. Is the code supposed to be this way?
Risposte (1)
Jan
il 27 Ago 2019
To accelerate your code, you can store the results of e.g. D1^5/pi^2 in a variable. This avoids 20 expensive power operations in the outer loop.
Matlab is optimized for vector and matrix operations. I assume the pile of different variables is less efficient.
Is there any hint, that the inner loop converges? In single typo can prevent eQ>tolQ from beeing reachable.
0 Commenti
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!