//Pitot Simulation //www.basicairdata.eu 2013 (c) JLJ clear ; //Setup variable graphics=0; //{1,0} graphics plots on/off //Pitot calibration factor c cpitot=0.9965790; //Model parameters R=1e-3; //[m] Internal radius of pressure tap tubing L=0.6;//[m] Line lenght Vu=50e-9;//[m^3] Sensor internal volume Vt=3.14*R^2*L;//[m^3] Line internal diameter sigma=0;// rhos=1.225; //[kg/m^3] Density of air at standard conditions mu=1.78938e-5//[kg/m/s] or [Pas]Aboslute fluid viscosity //mu=0; a0=340;//[m/s] speed of sound gammaratio=1.4;//Cp/Cv Specific heat ratio Pr=0.712;//Prandt Number nu=[1:0.5:1200]';//Frequency value range radiants [c r]=size(nu) for i=1:c alfa(i)=R*((rhos*nu(i)/mu)^0.5)*(complex(-2^0.5/2,2^0.5/2)) nj(i)=1.4; phaselag(i)=nu(i)/a0*(besselj(0,alfa(i))/besselj(2,alfa(i)))^0.5*((gammaratio/nj(i))^0.5) tfunction(i)=1/(cosh(phaselag(i)*L)+Vu/Vt*(sigma + 1/gammaratio)*nj(i)*phaselag(i)*L*sinh(phaselag(i)*L)) [phi(i),mag(i)]=phasemag(tfunction(i)) //mag(i)=10^(mag(i)/20) end fdt=frep2tf(nu/6.28,tfunction,2); fdt=factors(fdt) gain=syslin('c',1/0.9963631,1); //set the dc gain case 1 q=fdt*gain; //Pressure line transfer function if graphics==1 then for i=1:c [phi2(i),mag2(i)]=phasemag(repfreq(fdt,nu(i)/6.28)) end end //Launch xcos xcos('pitotdyn3.zcos') //Visualitation tasks if graphics==1 then xsetech([0,0,1,0.5]); plot2d(nu/6.28,mag',1) plot2d(nu/6.28,mag2',2) xtitle('Frequency response, single pressure tap line Original tf (black) Vs 2nd order identified tf (blue) www.basicairdata.eu JLJ') xlabel('Frequency [Hz]') ylabel('Magnitude OUT/IN [dB]') xgrid(1) xsetech([0.,0.5,1,0.5]); plot2d(nu/6.28,phi',1); plot2d(nu/6.28,phi2',2); xlabel('Frequency [Hz]') ylabel('Phase lag [Degree]') xgrid(1) end