+
+// This is a thermodynamically correct model of a mechanical vertical speed indicator:
+// It represents an aneroid in a closed (constant volume) casing with the aneroid internal pressure = static pressure
+// The casing has an orifice to static pressure
+// the mass flow through the orifice is calculated using compressible aerodynamics (but adiabatic and of course a perfect gas)
+// using the pressure in the casing and static pressure
+//
+// sadly at very low flows (small VS) in conjunction with the fact discrete timesteps (dt) are used, a numerical instability is formed.
+// this is counteracted by setting the massflow 0 at very small pressure differentials
+// this causes a small funny jump of your VSI when passing through 0...cannot be helped!
+//
+// also note the calibration is only valid for 0ft, so at higher altitudes, the vertical speed is not correct, but would indicate as a real mechanical VSI.
+// Only use for conventional mechanical VSI-s. Dont use in an Air Data Computer.
+//
+// (...and it is supposed to lag!)
+
+ _casing_airmass_kg = _casing_airmass_kg - _orifice_massflow_kgps * dt;
+ double new_density_kgpm3 = _casing_airmass_kg / Vol_casing;
+ _casing_pressure_Pa = _casing_pressure_Pa * pow(new_density_kgpm3 / _casing_density_kgpm3 , SG_gamma);
+ double casing_temperature_K = _casing_pressure_Pa / (new_density_kgpm3 * SG_R_m2_p_s2_p_K);
+
+ if( _casing_pressure_Pa - pressure_Pa > 0.0 ) {
+ Fsign = 1.0; //outflow, pos VS
+ } else {
+ Fsign = -1.0; //inflow, neg VS
+ }
+
+ if( fabs(_casing_pressure_Pa - pressure_Pa) < 0.01 ) {
+ orifice_mach = 0.0;
+ } else {
+ orifice_mach = sqrt(fabs (2.0*SG_cp_m2_p_s2_p_K / (SG_gamma * SG_R_m2_p_s2_p_K) * ( pow(pressure_Pa / _casing_pressure_Pa ,(SG_gamma-1)/SG_gamma ) -1 ) ) );
+ }
+
+ _orifice_massflow_kgps = Fsign * _casing_pressure_Pa / sqrt(casing_temperature_K) * sqrt(SG_gamma/SG_R_m2_p_s2_p_K) * orifice_mach * pow(1+(SG_gamma-1)/2*orifice_mach*orifice_mach,-(SG_gamma+1)/(2*(SG_gamma-1))) * A_orifice;
+
+ double vs_fpm = Fsign * sqrt( fabs( pressure_Pa - _casing_pressure_Pa ) ) * Factor_cal;
+ double vs_kts = vs_fpm / 60 * SG_FPS_TO_KT;
+ double vs_mps = vs_fpm / 60 * SG_FEET_TO_METER;
+
+ _speed_fpm_node
+ ->setDoubleValue(vs_fpm);
+ _speed_kts_node
+ ->setDoubleValue(vs_kts);
+ _speed_mps_node
+ ->setDoubleValue(vs_mps);
+
+ _casing_density_kgpm3 = new_density_kgpm3;
+