My code is running so long and never gives solution.
3 views (last 30 days)
Show older comments
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 Comments
Devineni Aslesha
on 26 Aug 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?
Answers (1)
Jan
on 27 Aug 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 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!