definite integration, I don't know if my integration is correct

1 次查看（过去 30 天）
Antonino Alesci
Antonino Alesci 2019-10-13
评论： Antonino Alesci 2019-10-13
I have this code, I have to integrate the expression expr, and I use polyfit and polyval to aquire the value of G that give to me the value T_Ninf = T. I think that my code is wrong, in particular the integration that I use is wrong because T_Ninf is a vector with high values.
clear all
close all
clc
%%% Temperature
T0 = 15; %°C
L = 0.0065; %°C/m (lapse rate)
h = 4000; %m
T_4000 = T0-(h*L);
T0k = 273.15+15;
T_4000k = T_4000+273.15; %K
%Pressione
P_4000_P0 = (T_4000k/T0k)^5.256;
P0 = 101325; %N/m^2
R = 287.3; %J/(kg K)
P_isa4000 = P_4000_P0*P0*exp((4000-0)/(287*T_4000k)); %N/m^2
%Densità
rho_4000 = P_isa4000/(R*T_4000k);
%Viscosità dinamica e cinematica
a = 1.46*10^-6; %coefficiente
b = 110; %K
mu = a*(T_4000k^1.5/(b+T_4000k));
nu = mu/rho_4000;
%Velocità del suono
gamma = 1.4;
a = sqrt(R*gamma*T_4000k);
%Stima di massima del numero di Reynold
R = (457*10^-3)/2; %m Il diametro
c = 0.2*R; %m valutazione al 75% del raggio
Omega = 561.82; %rad/s 5365 rpm
Re_075 = (rho_4000*Omega*0.75*R*0.2*R)/mu;
%Parametri
MM = 3.5; %kg;
M = MM/4;
% M_correct = M(8);
T = 9.81*M;
% T_correct = Tn(8);
A = R^2*pi; %area disco
r = linspace(10^-4,R,200); %valore estratto dal CAD
%Distribuzione radiale di circuitazione costante
% gam = (T*2*pi)/(A*rho_4000*Omega); %circuitazione
% u = gam./(2*pi*r); %velocità di rotazione al di sotto del disco
% u_inf = u/2;
% vi = sqrt(T/(2*rho_4000*A));
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Teoria dei vortici%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% numero infinito di pale %%%
G = linspace (10^-5,1,200);
r_s = r./(R.*G);
theta = acos((r_s.^6+3*r_s.^4+3*r_s.^2-1)./(r_s.^6+4*r_s.^4+3*r_s.^2+1));
omega_s = 6./(5+r_s.^2+2.*(1+r_s.^2).*cos(theta./3));
gam_opt = 2*pi*Omega.*omega_s.*r.^2; %distribusione ottima di circuitazione sul disco attuatore
u_r = gam_opt./(2*pi*r);
omega = u_r./r;
gam_opt_ad = omega_s.*(r./(R*G)).^2; %la velocità indotta è costante per un disco uniformemente caricato
% determinazione G di ottimo (ATTENZIONE !!!)
expr=@(r)2*pi*rho_4000.*(Omega-(omega./2)).*omega.*r.^3; %la funzione precedente è uguale, ma con omega esplicitato
T_Ninf = integral(expr,0,R,'ArrayValued',true);
p = polyfit(T_Ninf,G,2);
G_opt = polyval(p,T);
  2 个评论
Fabio Freschi
Fabio Freschi 2019-10-13
It seems correct: omega is an array. I tried to numerically calculate the integral, adding this piece of code at the bottom of yours
% sampling of r
rr = linspace(0,R,10000);
% loop over omega
for i = 1:length(omega)
expr1 = @(r)2*pi*rho_4000.*(Omega-(omega(i)./2)).*omega(i).*r.^3;
TT(i) = trapz(rr,expr1(rr));
end
In the end TT is practically identical to T_Ninf
Antonino Alesci
Antonino Alesci 2019-10-13
therefore, I think that somethink is wrong from the physical point of view. Anyway, thanks very much

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by