#define _USE_MATH_DEFINES
#include <math.h>
#include <stdio.h>
const double deg2rad= M_PI / 180;
double Avi2Trig(const double a) { return (180 - (a + 90)) * deg2rad; }
// --------------------------------------------------------------------------------------------------------------
const double dt= .001; // Passo d'integrazione
// Latitudine "media" per calcolo g
const double LAT= 39.7 * deg2rad, sl= sin(LAT), s2l= sin(2 * LAT);
const double aFronte= M_PI * 1*1 / 4., // diametro compressore= 1 m
m= 1800, // Aggiunti a 1500 kg di base la gondola, l'inversore e ~26 kg di fluidi (v. pagina "Caduta in mare")
// Durante la caduta del motore: 7,5e6 < Re < 16,7e6
// Re > 5e6 -> Cd= .52: G. Schewe: "On the force fluctuations acting on a circular cylinder in crossflow from subcritical up to transcritical Reynolds numbers"
// Re > 5e6 -> Cd= .52: NIEMANN & HOLSCHER, N. 1990 "A review of recent experiments on the flow past circular cylinders". Journal of Wind Engineering and Industrial Aerodynamics 33, 197 – 209
cAlato= .52 * (3.138 * 1.358); // Cr * A_laterale
// --------------------------------------------------------------------------------------------------------------
struct VECTOR { double x, y, z; };
double VMAG(VECTOR &v) { return sqrt(v.x*v.x + v.y*v.y + v.z*v.z); }
VECTOR _V(double x, double y, double z)
{
VECTOR vec = {x,y,z}; return vec;
}
VECTOR& operator+=(VECTOR& a,const VECTOR& b)
{
a.x += b.x; a.y += b.y; a.z += b.z;
return a;
}
VECTOR operator+ (const VECTOR &a, const VECTOR &b)
{
VECTOR c;
c.x = a.x+b.x; c.y = a.y+b.y; c.z = a.z+b.z;
return c;
}
VECTOR operator/ (const VECTOR &a, const double f)
{
VECTOR c;
c.x = a.x/f; c.y = a.y/f; c.z = a.z/f;
return c;
}
VECTOR operator* (const VECTOR &a, const double f)
{
VECTOR c;
c.x = a.x*f; c.y = a.y*f; c.z = a.z*f;
return c;
}
VECTOR operator* (const double f, const VECTOR &a)
{
return a*f;
}
// =========================================================
// === STATO ===
struct STS
{
VECTOR p, gs;
double tas, mach; // Informative
};
STS operator*(const double mul, const STS &s)
{
STS r; r.p= mul * s.p; r.gs= mul * s.gs;
return r;
}
STS operator+(const STS &s1, const STS &s2)
{
STS r; r.p= s1.p + s2.p; r.gs= s1.gs + s2.gs;
return r;
}
// ============================================================
// Assegna il meteo
void PARA(const double h, double &SKNT, double &DRCT, double &DEN, double &TEM)
{
const struct { double h, d, s, r, t; } // HGHT,DRCT,SKNT,DEN,TEMP
w[]= {
2000,319.0,13.0,1.1257,17.4,
3000,311.1,14.5,1.0924,16.0,
4000,302.4,15.8,1.0588,14.7,
5000,306.0,17.5,1.027,13.3,
6000,300.3,20.3,0.99593,12.0,
7000,294.1,23.5,0.96607,10.5,
8000,287.8,26.9,0.93673,8.9,
9000,281.7,30.4,0.90733,7.5,
10000,276.0,33.7,0.8784,6.1,
11000,272.6,37.6,0.85242,4.2,
12000,271.0,41.1,0.82701,2.4,
13000,269.3,44.8,0.80262,0.3,
14000,267.7,48.8,0.77853,-1.7,
15000,266.0,52.6,0.75437,-3.7,
16000,264.3,56.4,0.73017,-5.6,
17000,262.6,60.1,0.70621,-7.1,
18000,260.9,64.0,0.68425,-9.2,
19000,259.1,67.9,0.66335,-11.5,
20000,257.8,72.4,0.64142,-13.1,
21000,256.1,77.6,0.62047,-14.8,
22000,254.4,82.7,0.60027,-16.7,
23000,252.9,87.5,0.58091,-18.8,
24000,251.7,91.8,0.5623,-21.1,
25000,250.6,94.3,0.54471,-23.3,
26000,250.4,95.2,0.5279,-25.6
};
if(h < w[0].h) {
DRCT= w[0].d; SKNT= w[0].s; DEN= w[0].r; TEM= w[0].t;
return;
}
for(int i= sizeof(w) / sizeof(w[0]) - 1; i >= 0; i--) {
if(w[i].h < h) {
const double k= (h - w[i].h) / (w[i+1].h - w[i].h);
DRCT= w[i].d + (w[i+1].d - w[i].d) * k;
SKNT= w[i].s + (w[i+1].s - w[i].s) * k;
DEN = w[i].r + (w[i+1].r - w[i].r) * k;
TEM = w[i].t + (w[i+1].t - w[i].t) * k;
break;
}
}
}
// =================================================================================
STS DER(const STS &st)
{
double SKNT, DRCT, rho, TEM; PARA(st.p.z / .3048, SKNT, DRCT, rho, TEM);
// https://it.wikipedia.org/wiki/Accelerazione_di_gravit%C3%A0#Variazioni_locali_della_gravit%C3%A0_terrestre
double g= 9.7803184 * (1 + .0053024 * sl*sl - .0000059 * s2l*s2l) - 3.086e-6 * st.p.z;
double wv= SKNT * (1.852/3.6), wd= Avi2Trig(DRCT);
VECTOR w= _V(wv * cos(wd), wv * sin(wd), 0); // Vento [m/s]
VECTOR TAS= st.gs + w; double tas= VMAG(TAS);
// "Aircraft Performance Flight Testing" formula 4.30, pag. 57 - USATA dal 27-02-2023
double MACH= tas * 16.97498159 / (661.48*1.852/3.6 * sqrt(TEM + 273.15)); // 16.97... = sqrt(288.15)
/* MACH Cd+
.2 .175
.4 .18585
.8 .195
1 .196
*/
double Cadd;
if(MACH <= .4) Cadd= .175 + (.18585 - .175) / (.4 - .2) * (MACH - .2);
else if(MACH <= .8) Cadd= .18585 + (.195 - .18585) / (.8 - .4) * (MACH - .4);
else Cadd= .195 + (.196 - .195) / (1 - .8) * (MACH - .8);
double cFronte= .31 + Cadd; //printf("CB= %g kg/m^2, K= %g m^2/kg\n", m / cFronte / aFronte, cFronte * aFronte / m);
// Resistenza LATERALE DIVERSA da quella frontale
VECTOR F; // L'aereo tiene una prua = 180°
F.y= -.5 * rho * (st.gs.y + w.y) * tas * ( // Resistenza FRONTALE del motore
cFronte * aFronte // Windmilling
+
.0025 * M_PI * 5.255 * (1.143 + 1.524 + 0.746) / 3 // Resistenza viscosa ("Gas Turbine Performance" F5.5.1 e F5.5.2)
);
F.z= -.5 * rho * st.gs.z * tas * (
cFronte * aFronte
+
.0025 * M_PI * 5.255 * (1.143 + 1.524 + 0.746) / 3
);
// Resistenza LATERALE del motore
F.x= -.5 * rho * (st.gs.x + w.x) * tas * cAlato;
F.z -= m * g; // Gravità
STS der;
der.p= st.gs; // Derivata della posizione: p' = v
der.gs= F / m; // Derivata della velocità: v' = F / m
der.tas= tas; der.mach= MACH; // Informative
return der;
}
// =================================================================================================================
int main(void)
{
// Stato iniziale del motore
const double gs0= 466.98687690461679 * 1.852/3.6,
rotta= Avi2Trig(167.79960430883506);
STS st= {
0, 0, 25926.956176757812 * .3048, // Posizione [m]
gs0 * cos(rotta), gs0 * sin(rotta), 0 // Velocità inerziale [m/s]
};
int c_stp= 0; // Numero di passi
double h_old= HUGE_VAL, h0= st.p.z, t;
for(t= 0; ; t += dt) {
// RK-2
STS f1 = DER(st);
STS der= DER(st + .5 * dt * f1);
c_stp++;
// Mostra l'integrazione ogni 50 m o quando la quota è inferiore a 1 m
if(h_old - st.p.z > 50 || st.p.z < 1) {
printf("%5d) %8.4f, h= %7.1f, P: %6.1f %7.1f %7.1f, GS: %4.1f %6.1f %6.1f -> %6.1f, a: %5.2f %5.2f %5.2f -> %5.2f",
c_stp, t, st.p.z / .3048,
st.p.x, st.p.y, st.p.z, // POS
st.gs.x, st.gs.y, st.gs.z, VMAG(st.gs), // VEL
der.gs.x, der.gs.y, der.gs.z, VMAG(der.gs) // ACC
);
printf(", TAS= %.1f kn %.3f mach\n", der.tas * 3.6/1.852, der.mach);
h_old= st.p.z;
}
// Motore arrivato al suolo, interrompe l'integrazione
if(st.p.z < 0) break;
// Fa avanzare lo stato del motore
st = st + dt * der;
}
printf("\nt: %f s\n", t);
printf("POS [m ]: x= %f, y= %f, z= %f\n", st.p.x, st.p.y, st.p.z);
printf("VEL [m/s]: x= %f, y= %f, z= %f\n", st.gs.x, st.gs.y, st.gs.z);
double dst= sqrt(st.p.x*st.p.x + st.p.y*st.p.y);
printf("\nORTODROMIA da POS iniziale a finale: d= %f nm (%.1f m), al= %f deg\n", dst / 1852, dst, 90 - atan2(st.p.y, st.p.x) / deg2rad);
return 0;
}
