Generate Random Numbers to satisfy a variable on the left and right hand side
2 次查看(过去 30 天)
显示 更早的评论
Greetings and Hi !
Supposedly if i would like to generate a random number for M1 and M2 which would then be inputted into the coming equations.
R = 287; %Gas constant
T01= 339; %Stagnation Temperature or total inlet temperature
Patm = 101325;
P01 = [124638.717 126432.080 128225.442 130018.805 131812.168 135398.894 137192.256637168 142572.345132743 142572.345132743 147952.433628319 155125.884955752 165886.061946903 167679.424778761 183819.690265487 201753.318584071 225067.035398230 234033.849557522 239413.938053097 248380.752212389 253760.840707965 263445 273577.500000000 283710 293842.500000000 303975 309041.250000000 310054.500000000 310561.125000000 311067.750000000 311574.375000000 314107.500000000];
A1 = 0.0025686; %Area inlet volute
kc = 1.4; %specific heat ratio
kc0 = 0.2; %specific heat ratio at total state
kc1 = 2.498; %specific heat ratio at inlet volute
kc2 = 3.498; %specific heat ratio at inlet rotor
%Mach Number 1 (Guess)
M1 = rand(1,31);% enter the vector M1
rho01 = P01./(R*T01);
rho1 = rho01./(1 + kc0 * M1.^2).^kc1;
for iter_1 = 1:length(M1)
rho01 = P01./(R*T01); %Stagnation density
rho1 = rho01./(1 + kc0 * M1.^2).^kc1; %flow density inlet volute
T1 = T01./(1 + kc0 * M1.^2); %static temperature inlet volute
P1 = P01./(1 + kc0 *M1.^2).^kc2; %static pressure inlet volute
C1 = M1.*(kc*R*T1).^0.5; %absolute velocity inlet
mdot1 = rho1.*A1.*C1; %mass flow rate inlet volute
end
%% Station 2
M2 = rand(1,31);% enter the value M2
%M2 = [0.320 0.333 0.346 0.359 0.371 0.395 0.407 0.441 0.441 0.473 0.514 0.571 0.580 0.660 0.744 0.848 0.888 0.912 0.953 0.977 1.023 1.074 1.130 1.192 1.271 1.324 1.336 1.344 1.351 1.359 1.411]
T02 = 339; %Stagnation temperature or total temperature at rotor inlet
N = 47503; %Rotor speed
omega = (2*pi/60)*N; %Rotor 100% speed angular velocity
ktp_volute= 0.9; %Coeeficient Total Pressure Loss Volute
SW = 0.95; % Swirl Coefficient
B2 = 1; %Inlet rotor blockage factor
r1 = 0.0739; %Inlet volute radius
r2max = 0.047576; %Rotor shroud tip radius at leading edge
r2min = 0.036006; %Rotor hub radius at leading edge
rrms2 = (r2max + r2min)/2; %Average between r2max and r2min
U2 = rrms2 * omega; %Blade speed rotor tip
for iter_2 = 1:length(M2)
Ctheta2_exitvolute = C1.*SW.*(r1 / rrms2 ); %Tangential absolute velocity inlet rotor
P2 = P01./(((1 + kc0 * M2.^2).^kc2).*(1+ktp_volute)-ktp_volute); %Static pressure inlet rotor
T2 = T02./(1 + kc0 * M2.^2); %static temperature inlet rotor
P02 = P2.*(1 + kc0 * M2.^2).^kc2; %Total pressure inlet rotor
rho2 = P2./(R*T2); %Stagnation density
% using continuity equation, mass flow rate inlet rotor = mass flow rate inlet volute
% mdot1 = mdot2
mdot2 = mdot1; %mass flow inlet rotor
leh = 0.018; %rotor leading height
A2 = pi * (r2max +r2min) * leh;%Area Inlet Rotor
%Absolute velocity inlet rotor
alpha2 = (atand((rho2./rho1).*((A2/rrms2)./( A1/r1)).*SW));
C2 = M2.*(kc*R*T2).^0.5; %absolute velocity inlet volute
Rad = deg2rad(alpha2);
Ctheta2_inletrotor = C2.*sin(Rad); %Tangential absolute velocity inlet rotor
Ctheta2_exitvolute = Ctheta2_inletrotor ; %DeltaCtheta2
end
However what i need from M1 and M2 is to generate numbers in which when inputted will give the same value between Ctheta2_exitvolute and Ctheta2_inletrotor
5 个评论
Torsten
2022-3-24
Should the relation between M1 and M2 to satisfy
Ctheta2_exitvolute = Ctheta2_inletrotor
depend on P01 ?
回答(2 个)
Torsten
2022-3-25
编辑:Torsten
2022-3-25
Here is a code with which you get M2 versus M1 curves for the vector P01 of pressures for which the relation
Ctheta_exitvolute - Ctheta_inletrotor = 0
is satisfied.
Interestingly, the curves don't depend on P01. That's why I asked whether the M1-M2 relation should depend on P01.
Be a little patient - the code takes quite a while.
P01 = [124638.717 126432.080 128225.442 130018.805 131812.168 135398.894 137192.256637168 142572.345132743 142572.345132743 147952.433628319 155125.884955752 165886.061946903 167679.424778761 183819.690265487 201753.318584071 225067.035398230 234033.849557522 239413.938053097 248380.752212389 253760.840707965 263445 273577.500000000 283710 293842.500000000 303975 309041.250000000 310054.500000000 310561.125000000 311067.750000000 311574.375000000 314107.500000000];
m1 = 0:0.001:0.6;
m2 = 0:0.001:5;
for k=1:numel(P01)
p01 = P01(k);
for i=1:numel(m1)
M1 = m1(i);
resl = fun(m2(1),M1,p01);
flag = 0;
for j=1:numel(m2)-1
resr = fun(m2(j+1),M1,p01);
if resl*resr <=0
sol_m2(i,k) = (m2(j)+m2(j+1))/2;
flag=1;
break
end
resl = resr;
end
if flag==0
sol_m2(i,k)=NaN;
end
%if mod(i,100)==0
% i
%end
end
k
end
plot(m1,sol_m2)
end
function res = fun(x,M1,P01)
R = 287; %Gas constant
T01= 339; %Stagnation Temperature or total inlet temperature
Patm = 101325;
%P01 = [124638.717 126432.080 128225.442 130018.805 131812.168 135398.894 137192.256637168 142572.345132743 142572.345132743 147952.433628319 155125.884955752 165886.061946903 167679.424778761 183819.690265487 201753.318584071 225067.035398230 234033.849557522 239413.938053097 248380.752212389 253760.840707965 263445 273577.500000000 283710 293842.500000000 303975 309041.250000000 310054.500000000 310561.125000000 311067.750000000 311574.375000000 314107.500000000];
A1 = 0.0025686; %Area inlet volute
kc = 1.4; %specific heat ratio
kc0 = 0.2; %specific heat ratio at total state
kc1 = 2.498; %specific heat ratio at inlet volute
kc2 = 3.498; %specific heat ratio at inlet rotor
rho01 = P01./(R*T01);
rho1 = rho01./(1 + kc0 * M1.^2).^kc1;
T1 = T01./(1 + kc0 * M1.^2); %static temperature inlet volute
P1 = P01./(1 + kc0 *M1.^2).^kc2; %static pressure inlet volute
C1 = M1.*(kc*R*T1).^0.5; %absolute velocity inlet
mdot1 = rho1.*A1.*C1; %mass flow rate inlet volute
T02 = 339; %Stagnation temperature or total temperature at rotor inlet
N = 47503; %Rotor speed
omega = (2*pi/60)*N; %Rotor 100% speed angular velocity
ktp_volute= 0.9; %Coeeficient Total Pressure Loss Volute
SW = 0.95; % Swirl Coefficient
B2 = 1; %Inlet rotor blockage factor
r1 = 0.0739; %Inlet volute radius
r2max = 0.047576; %Rotor shroud tip radius at leading edge
r2min = 0.036006; %Rotor hub radius at leading edge
rrms2 = (r2max + r2min)/2; %Average between r2max and r2min
U2 = rrms2 * omega; %Blade speed rotor tip
Ctheta2_exitvolute = C1.*SW.*(r1 / rrms2 ); %Tangential absolute velocity inlet rotor
P2 = P01./(((1 + kc0 * x.^2).^kc2).*(1+ktp_volute)-ktp_volute); %Static pressure inlet rotor
T2 = T02./(1 + kc0 * x.^2); %static temperature inlet rotor
P02 = P2.*(1 + kc0 * x.^2).^kc2; %Total pressure inlet rotor
rho2 = P2./(R*T2); %Stagnation density
% using continuity equation, mass flow rate inlet rotor = mass flow rate inlet volute
% mdot1 = mdot2
mdot2 = mdot1; %mass flow inlet rotor
leh = 0.018; %rotor leading height
A2 = pi * (r2max +r2min) * leh;%Area Inlet Rotor
%Absolute velocity inlet rotor
alpha2 = (atand((rho2./rho1).*((A2/rrms2)./( A1/r1)).*SW));
C2 = x.*(kc*R*T2).^0.5; %absolute velocity inlet volute
Rad = deg2rad(alpha2);
Ctheta2_inletrotor = C2.*sin(Rad); %Tangential absolute velocity inlet rotor
res = Ctheta2_exitvolute - Ctheta2_inletrotor;
end
4 个评论
Walter Roberson
2022-4-20
However, I am still unable to achieve the desired value 0
It is possible that your goal cannot be reached. You are using floating point expressions, and calculating values in two different ways, but you are expecting the two values to be bit for bit identical . That is what it means to compare A == B -- that A and B are bit for bit identical . Exact equality.
For example, adjacent representable numbers
format long g
A = hex2num('3fe63c093c802fe9')
B = hex2num('3fe63c093c802fe8')
fprintf('A = %.15g\n', A)
fprintf('B = %.15g\n', B)
They look the same
A == B
But they compare different
A-B
There is an actual difference
fprintf('%.99g\n', A)
fprintf('%.99g\n', B)
Because of this, you should never be using == on floating point numbers that are calcuated in different ways. Instead you should be calculating
if abs(A - B) < SomeTolerance
Torsten
2022-4-20
编辑:Torsten
2022-4-20
ok.
This code should return M1, given M2 and a starting guess for M1, such that the pair (M1,M2) satisfies
Ctheta2_exitvolute == Ctheta2_inletrotor
M2 = 2.0;
M10 = 0.5;
M1 = fsolve(@(x)fun(x,M2),M10)
fun(M1,M2)
function res = fun(m1,M2)
%% Calculation
%Station 1
R = 287; %Gas constant
T01= 339; %Stagnation Temperature or total inlet temperature
Patm = 101325;
P01 = 124638.717;
A1 = 0.0025686; %Area inlet volute
kc = 1.4; %specific heat ratio
kc0 = 0.2; %specific heat ratio at total state
kc1 = 2.498; %specific heat ratio at inlet volute
kc2 = 3.498; %specific heat ratio at inlet rotor
%% Station 2
%M2 = rand(1)
T02 = 339; %Stagnation temperature or total temperature at rotor inlet
N = 47503; %Rotor speed
omega = (2*pi/60)*N; %Rotor 100% speed angular velocity
ktp_volute= 0.9; %Coeeficient Total Pressure Loss Volute
SW = 0.95; % Swirl Coefficient
B2 = 1; %Inlet rotor blockage factor
r1 = 0.0739; %Inlet volute radius
r2max = 0.047576; %Rotor shroud tip radius at leading edge
r2min = 0.036006; %Rotor hub radius at leading edge
rrms2 = (r2max + r2min)/2; %Average between r2max and r2min
U2 = rrms2 * omega; %Blade speed rotor tip
%for iter_2 = 1:length(M2)
P2 = P01./(((1 + kc0 * M2.^2).^kc2).*(1+ktp_volute)-ktp_volute); %Static pressure inlet rotor
T2 = T02./(1 + kc0 * M2.^2); %static temperature inlet rotor
P02 = P2.*(1 + kc0 * M2.^2).^kc2; %Total pressure inlet rotor
rho2 = P2./(R*T2); %Stagnation density
leh = 0.018; %rotor leading height
A2 = pi * (r2max +r2min) * leh;%Area Inlet Rotor
%m1 = rand() %Changing M1 to get value DeltaCtheta2
%for iter_1 = 1:length(m1)
rho01 = P01./(R*T01); %Stagnation density
rho1 = rho01./(1 + kc0 * m1.^2).^kc1; %flow density inlet volute
T1 = T01./(1 + kc0 * m1.^2); %static temperature inlet volute
P1 = P01./(1 + kc0 *m1.^2).^kc2; %static pressure inlet volute
C1 = m1.*(kc*R*T1).^0.5; %absolute velocity inlet
mdot1 = rho1.*A1.*C1; %mass flow rate inlet volute
Ctheta2_exitvolute = C1.*SW.*(r1 / rrms2 ); %Tangential absolute velocity inlet rotor
% using continuity equation, mass flow rate inlet rotor = mass flow rate inlet volute
% mdot1 = mdot2
mdot2 = mdot1; %mass flow inlet rotor
%Absolute velocity inlet rotor
alpha2 = (atand((rho2./rho1).*((A2/rrms2)./( A1/r1)).*SW));
Rad = deg2rad(alpha2);
%Absolute velocity inlet rotor
%C2 = M2(iter_2).*(kc*R*T2).^0.5;
C2 = M2.*(kc*R*T2).^0.5;
Ctheta2_inletrotor = C2.*sin(Rad);
res = Ctheta2_exitvolute - Ctheta2_inletrotor;
%DeltaCtheta
%syms Ctheta2_exitvolute Ctheta2_inletrotor
%DeltaCtheta2 = solve([Ctheta2_exitvolute - Ctheta2_inletrotor == 0],Ctheta2_exitvolute)
%tolerance = Ctheta2_exitvolute - Ctheta2_inletrotor
%while isempty(tolerance) || tolerance == 0 || tolerance > 1
%C1 = m1.*(kc*R*T1).^0.5 ;
%Ctheta2_exitvolute = (m1.*(kc*R*T1).^0.5).*SW.*(r1 / rrms2 ) ;
%accuracy = m1; %abs(Ctheta2_exitvolute - Ctheta2_inletrotor)
%while tolerance ~= 0
%m1 = m1 + 0.01; %want to continue iterate until achieve a value that corresponds with DeltaCtheta2
%C1 = m1.*(kc*R*T1).^0.5;
%accuracy = abs(m1);
%if (Ctheta2_exitvolute < Ctheta2_inletrotor)
%m1 = (m1 + (m1+0.01))/2
%Ctheta2_exitvolute = m1.*(kc*R*T1).^0.5
%fprintf('Criteria NOT satisfied with value Mach Number 1 of %f\n', m1);
%elseif Ctheta2_exitvolute > Ctheta2_inletrotor
%m1 = (m1 + (m1 - 0.01))/2
%Ctheta2_exitvolute = m1.*(kc*R*T1).^0.5;
%fprintf('Criteria NOT satisfied with value Mach Number 1 of %f\n', m1);
%else abs(Ctheta2_exitvolute == Ctheta2_inletrotor) %<= eps(Cheta2_inletrotor)
%Continue
%fprintf('Criteria satisfied with value Mach Number 1 of %f\n', m1);
%break
%end
%C1
%m1
%tolerance
%end
%if Ctheta2_exitvolute == Ctheta2_inletrotor
%fprintf('Delta_Ctheta2 satisfied with value of %f\n', Ctheta2_exitvolute, Ctheta2_inletrotor);
% end
%end
%end
%end
end
8 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 PID Controller Tuning 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!