-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmain.m
155 lines (153 loc) · 6.9 KB
/
main.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
%%% Pipe optimal design: smaller diameter / smaller pressure drop
% CONTENT
% This script solves the problem of finding the optimal diameter for a
% straight pipe coneying an incompressible fluid, for different mass flow
% rates.
% The best design solution is considered the one that minimizes both
% the pressure drop per unitary pipe length and the pipe diameter.
%
% DESIGN VARIABLE
% There is 1 design variable (the pipe diameter) and 2
% competitive objectives, so the problem is solved with the Pareto-optimum
% approach, that will output a compromise solution.
%
% The user can express both lower and upper boundaries for the search of
% the optimal pipe diameter and can provide lower and upper boundary for
% the flow velocity [m/s], thus constraining the optimization problem.
%
% PROBLEM CONSTANTS
% To perform the optimization, some fluid properties such as fluid dynamic
% viscosity and density shall be imposed. The rated mass flow rate to be
% conveyed by the pipe and the pipe absolute roughness (that depends from
% the pipe material and conditions) shall be given as constant inputs.
%
% SOLVER CHOICE
% Non dominated configurations are found with the multi-objective genetic
% algorithm solver "gamultiobj", but other algorithms such as
% "paretosearch" can be used, editing the solver section.
%
% Note that both multi-objective solvers do not support inter variables as
% default, while "ga" solver enables also mixed-integer problems, but was
% conceived for mono-objective optimization only. One could also define a
% scalar overall cost function, like [1], but the purpose of this study is
% to consider the 2 objectives strictly separated.
%
% HOW IT WORKS
% For each mass flow rate, the program calculates non-dominated designs in
% terms of pipe diameter [m] and pressure drop per meter [mm H20/m].
% The best solution is chosen with the TOPSIS method.
% For a closer look on this outranking method, a good explanation is
% provided by [2]. Since there are only 2 objectives, I chose to assignan
% equal weight to each one, thinking that they have the same importance.
% The user is free to change the weights when calling the topsis.m
% function. For example, one could assign w=[0.25, 0.75] to make the
% diameter objective more relevant.
%
% Finally, the best design for each flow rate is represented in the "pipe
% chart", that the log-log chart commonly used by technicians to speed up
% pipelines sizing.
%
% CHANGE THE FLUID TYPE
% This is an interesting direction for investigation. One should remember
% that the practical log-log charts for pipe sizing are developed for a
% single type of pipe (i.e. carbon steel or plastic pipes), expressing a
% reasonable value of epsilon and operated with a well defined working
% fluid (i.e. water or air are the most common). In few words, these charts
% are not universal, but they are drawn depending from
% 1. pipe material
% 2. fluid properties (density & viscosity)
%
% Try to input the values of a tipical lubrication oil at design
% temperature of 40 °C, which are:
% - kinematic viscosity @ 40: 140 cSt
% - density @ 40: 950 kg/m^3
% that give dynamic visocity of 144*0.95*10^-6 [Poiseuille = Pa*s]
% use a plastic polymer pipe, with epsilon = 0.002*1e-3 [m] and feasable
% oil speed between 0.05 and 2 [m/s].
%
% Or try to solve the same problem for a fuel oil. For example consider:
% - kinematic viscosity @80: 985 kg/m^3
% - density @ 80: 985 cSt
% - high quality steel pipe, epsilon=0.01*1e-3;
% - admitted fluid velocity between 0.01 and 1.5 m/s
%
% References:
%
% [1] S. B. Genic, B. M. Jacimovic, V. B. Genic, Economic optimiztion of
% pipe diameter for complete turbulence, Energy and Buildings
% 45 (2012) p. 335–338. doi:10.1016/j.enbuild.2011.10.054
%
% [2] X. Li et al. Application of the Entropy Weight and TOPSIS Method in
% Safety Evaluation of Coal Mines, Procedia Engineering.
% doi:10.1016/j.proeng.2011.11.2410
%
clearvars;
clc;
%% input data
rho=1000; % fluid density constant [kg/m3] e.g. for water @ 20 °C
mu=9*10^-4; % fluid dynamic viscosity [Pa*s]
% this can also be expressed as mu=nu*rho, if the dynamic
% viscosity is given. pay attention that these values
% depend from the fluid temperature.
roughness=0.2*1e-3; % pipe rated absolute roughness [m]
% 0.001*1e-3 for floor heating pipe [m]
% 0.05*1e9-3 steel pipe with fine surface [m]
% 0.5*1e-3 common steel pipe [m]
%% define variable boundaries
D_min=5*1e-3; % lower bound for the pipe diameter in [m]
D_max=1000*1e-3; % upper bound for the pipe diameter in [m]
%% define constraints on the fluid velocity
% max & min fluid velocity allowed [m/s].
% These are common values that can be retrieved from technical manuals.
% The user can adjust as required.
velocity_boundaries=[0.1, 3]; % velocity [m/s]
%% input mass flow rates array
% Estimate best design for all these flow rates.
% m_dot_plot=[0.2, 0.3, 0.4, 0.5, 0.5, 0.7, 0.8, 0.9, 1, 2, 3, 4, 5, ...
% 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300, 400];
m_dot_plot=linspace(0.1, 1, 10);
%% common diameters array
% This array is needed to plot the pipes chart. No influence on the
% optimization problem.
D_plot= [15, 20, 30, 40, 50, 65, 80, 100, ...
125, 150, 200, 250, 300, 350, 400]./1000;
%% solve the problem
all_results=table();
for k=1:numel(m_dot_plot)
m_dot=m_dot_plot(k);
[results{k}, exitflag(k)]=optimize_pressure_drop_with_pipe_diameter(m_dot, ...
rho, mu, roughness, D_min, D_max, velocity_boundaries);
% these lines to give a feedback on how the calculation is proceeding.
if exitflag(k) > 0
message=sprintf('run %d of %d successful',k, numel(m_dot_plot));
% sum up all non-dominated designs in a cell arra, one for each
% mass flow rate
all_results(k,:)=results{k}(1,:);
else
message=sprintf('run %d of %d failed',k, numel(m_dot_plot));
end
disp(message);
end
%% save the best designs for each flow rate
writetable(all_results, 'optimal_pipe_design.csv');
%% plot the all optimal design in the "pipes chart"
figure(1)
dp_name='dp [mm H2O/m]';
d_name='d [mm]';
% plot the "pipes chart"
if numel(m_dot_plot) > 1
plot_pipes_chart(m_dot_plot, rho, mu, D_plot, roughness, ...
[0.01, 1000], [0.1, max(m_dot_plot)]);
%
% plot the best designs
for k=1:numel(m_dot_plot)
loglog(all_results.(dp_name)(k), all_results.m_dot(k),'xk');
hold on;
text(all_results.(dp_name)(k), all_results.m_dot(k), ...
num2str(all_results.(d_name)(k)),'VerticalAlignment','bottom', ...
'HorizontalAlignment','left');
end
% savefig('pipes_chart_with_optimal_designs.fig');
else
disp('pipes chart will not be plotted for only 1 flow rate point');
end