-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathbag_flow_plot.m
144 lines (116 loc) · 3.56 KB
/
bag_flow_plot.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
%% Flow calculations
%constants of the system:
%environmental constants
pa = 101325; %atmpsoheric pressure [kPa]
gamma = 1.4; %specific heat ratio [1]
T = 300; %outside temperature [K]
R = 287; %air gas constant [J/(kg*K)]
g = 9.81; %gravitational acceleration m/s^2
N = 4; %number of skates
%Pod 2
%
% %pod/skate properties
% M = 909; %mass of pod, [kg]
% %M = 324.319;
%
% %length/width/area
% L = 0.9144; % length of skate [m]
% %W = 0.127; % width of skate [m]
% W = 0.3048; %skate width [m]
M = 5;
L = 0.1397;
W = 0.0762;
A_skate = W*L; %skate area, [m^2]
Per_skate = 2*(W+L); %perimeter, [m]
porosity = 1;
%porosity = 0.04;
A = porosity * A_skate*N;
Per = Per_skate*N;
h_skate = 0.0445; %bezel height
h_bag = 0.020066;
%h = h_skate + h_bag;
h = h_bag;
V = W*L*h;
mu = 1.8e-5; %dynamic viscosity [Pa * s]
%Find equilibrium states
%each skate needs to hold up 1/4 of pod's weight
zstar = 5e-5;
%input equilibrium
%mass flow at equilibrium
min0 = obj.porosity * obj.A * obj.M * obj.g / (obj.mu * obj.h_bag);
%find bottom pressure from mass flow
mass_escape_balance = @(p) mass_flow_out(p, obj.pa, z0, obj.gamma, obj.Per, obj.R, obj.T) - min0;
p0 = fzero(mass_escape_balance, obj.pa);
%bag pressure to keep the skate aloft
pb0 = p0 + obj.M*obj.g/obj.A;
pstar = pa + (M*g)/A; %skate pressure, Pa
minstar = porosity * m_out(pstar, pa, zstar, gamma, Per, R, T);
%start mass flow calculations over grid
Ng = 71;
%z = linspace(0, 0.005, N);
z = linspace(0, 5*zstar, Ng);
p = linspace(pa, pstar + 2.5*(pstar - pa), Ng);
[zz, pp] = meshgrid(z, p);
mass_out = m_out( pp, pa, zz, gamma, Per, R, T );
%unit conversions
kgs_to_scfm = 1800.24; %kg/s to standard cubic feet per minute
flow_out = mass_out .* kgs_to_scfm;
fstar = minstar*kgs_to_scfm;
%% Visualizations
%contour plot 2d + 3d surface
figure(48)
clf
%surf(pp, zz, mass_out)
hold on
%contour(pp, zz, flow_out, 15, 'ShowText', 'on');
surf(pp, zz, flow_out, 'EdgeColor', 'None');
alpha 0.5
contour(pp, zz, flow_out, 20);
plot3([pa, pstar], [zstar, zstar], [fstar, fstar], '--k')
plot3([pstar, pstar], [0, zstar], [fstar, fstar], '--k')
plot3([pstar, pstar], [zstar, zstar], [0, fstar], '--k')
scatter3(pstar, zstar, fstar, 300, '*r')
hold off
xlabel('skate pressure (Pa)')
ylabel('ride height (m)')
zlabel('skate volume flow out (scfm)')
title('Q/skate vs. z and p')
c = colorbar;
c.Label.String = 'Per Skate Flow (scfm)';
view(3)
%contour plot 2d only
figure(49)
clf
hold on
%surf(pp, zz, flow_out)
contour(pp, zz, flow_out, 10, 'ShowText', 'on')
scatter(pstar, zstar, 100, '*r')
txtZ = strcat('z^* = ', num2str(zstar*1000, 2), ' mm');
txtP = strcat('p^* = ', num2str(pstar/1000, 4), ' kPa');
txtQ = strcat('Q_{in} = ', num2str(fstar, 4), ' scfm');
text(pa, zstar - 3e-5, txtZ)
text(pstar + 1e2, 0+ 3e-5, txtP)
text(pstar + 2e2, zstar,txtQ)
plot([pstar, pstar], [0, zstar], '--k')
plot([pa, pstar], [zstar, zstar], '--k')
hold off
xlabel('skate pressure (Pa)')
ylabel('ride height (m)')
zlabel('skate volume flow out (scfm)')
title('Q/skate vs. z and p')
c = colorbar;
c.Label.String = 'Per Skate Flow (scfm)';
axis square
function [ mass_out ] = m_out( pp, pa, zz, gamma, Per, R, T )
%m_out finds the mass flow rate out in an air skate in kg/s
%pp: skate pressure
%pa: atmospheric (outside) pressure
%zz: current ride height
%gamma: specific heat ratio, ~1.4
%Per: total perimeter of air skates
%R: air gas constant
%T: atmospheric temperature
outside_term = pp.*Per.*zz/sqrt(R*T);
vout2_term = 2*gamma/(gamma-1) * ((pa./pp).^(2/gamma) - (pa./pp).^((gamma+1)/gamma));
mass_out = outside_term .* sqrt(vout2_term);
end