Spadanie swobodne z uwzględnieniem oporu powietrza.

Do napisania tego tekstu zdopingowała mnie dyskusja na pl.rec.paralotnie, pierwotnie o opadaniu na zapasie ... i zaproponowane przez jednego z kolegów doświadczenie, polegające na jednoczesnym upuszczeniu w pozycji pionowej, dwóch identycznych plastikowych butelek:

- pierwszej z zawartością 0,2 l wody
- drugiej z 1,3 litra wody.

Wyobrażenie sobie przebiegu zjawiska fizycznego w czasie może być trudne, nawet jak zna się opisujące je wzory fizyczne. Dlatego napisałem prosty kawałek kodu modelujacego spadanie 2 ciał o różnych masach ale jednakowym kształcie. W kolejnych krokach obliczeń siła działajaca na ciało (i powodująca jego przyspieszanie) to jego ciężar pomniejszony o opór powietrza. Opór powietrza jest liczony na podstawie prędkości z poprzedniego kroku, niedokładności wynikające z tego uproszczenia można zmniejszać zmniejszając krok obliczeniowy. Do celów demonstracji wystarczająco dobre wyniki moża uzyskać przy kroku 0,5-0,1 sekundy. Wydaje mi się, że ten poziom matematycznego opisu i jego implementacja w taki prosty sposób powinny być zrozumiałe dla każdego ...

Przyjęto pole poprzecznego przekroju 0,006 metra kwadratowego (pole koła o średnicy 9 cm), stały współczynnik oporu 0,1 niezależnie od prędkości (wzrastająca prędkość w wpływa na skalę zjawiska i opisującą ją liczbę Re ale jak to wpływa na opór butelki ?) i stałą gęstość powietrza (gęstość rośnie wraz ze spadkiem wysokości), czas eksperymentu 10 sekund.

Kod można uruchomić w internetowej 'piaskownicy programistycznej' takiej jak ideone.com (ten kod pod linkem https://ideone.com/CCEYnT ) lub w dowolnym innym środowisku z językiem C, a modyfikując odpowiednio parametry (np. masy ciał, pole poprzecznego przekroju, współczynnik oporu, ...) wirtualnie przeprowadzać różne wersje doświadczenia, wyniki w postaci tabelki ( rysowanie wykresów w wersji 2.0, animacje w 3.5 ... ).

#include < stdio.h >

// stałe fizyczne
const double ro = 1.25; // gestosc powietrza [kg/m^3]
const double g  = 9.81; // przyspieszenie ziemskie [m/s^2]

double Faero(double S, double C, double v) { return v*v*S*C*ro/2; }

int main(void) {
// parametry eksperymentu stale w czasie
double m1 = 0.2; // masa ciala 1 [kg]
double m2 = 1.3; // masa ciala 2 [kg]

double Cx1  =  0.1; // wspolczynnik oporu ciala 1
double Cx2 =   0.1;

double S_1 = 0.006; // pole poprzecznego przekroju [m^2]
double S_2 = 0.006;

double dt    = 0.1;  // modelowy odcinek czasu [s]
double Tmax = 10.0;  // limit czasu [s]

// zmienne w czasie
double a1 = 0.0; // przyspieszenie [m/s^2]
double a2 = 0.0;
double f1 = 0.0; // sila oporu [N]
double f2 = 0.0;
double v1 = 0.0; // predkosc [m/s]
double v2 = 0.0;
double s1 = 0.0; // droga [m]
double s2 = 0.0;

for(double t=0; t < Tmax; t+=dt) {
f1 = Faero(S_1, Cx1, v1); // opor liczony dla predkosci z poprzedniego kroku
f2 = Faero(S_2, Cx2, v2);
a1 = (m1*g - f1) / m1; // przyspieszenie = (ciezar - opor) / mase
a2 = (m2*g - f2) / m2;
double v1p = v1; // predkosc poczatkowa w kroku - do policzenia sredniej
double v2p = v2;
v1 += a1*dt; // predkosc wzrasta o 'przyspieszenie * czas_kroku'
v2 += a2*dt;
s1 += (v1p + v1)/2 * dt; // droga wzrasta o 'sredna_predkosc * czas_kroku'
s2 += (v2p + v2)/2 * dt;
printf("%3.2fs f[%6.2f,%6.2f] a[%4.2f,%4.2f] v[%5.2f,%5.2f] s[%5.2f,%5.2f]\n",
       t+dt, f1,f2, a1,a2, v1,v2, s1,s2);
}

return 0;
}

Wyniki [pierwsze_ciało, drugie_ciało]:

0.50s f[  0.01,  0.01] a[9.78,9.81] v[ 4.90, 4.90] s[ 1.23, 1.23]
0.60s f[  0.01,  0.01] a[9.76,9.80] v[ 5.88, 5.88] s[ 1.76, 1.77]
0.70s f[  0.01,  0.01] a[9.75,9.80] v[ 6.85, 6.86] s[ 2.40, 2.40]
0.80s f[  0.02,  0.02] a[9.72,9.80] v[ 7.82, 7.84] s[ 3.13, 3.14]
0.90s f[  0.02,  0.02] a[9.70,9.79] v[ 8.79, 8.82] s[ 3.97, 3.97]

2.90s f[  0.26,  0.28] a[8.52,9.60] v[27.13,28.24] s[40.29,41.10]
3.00s f[  0.28,  0.30] a[8.43,9.58] v[27.97,29.19] s[43.05,43.97]
3.10s f[  0.29,  0.32] a[8.34,9.56] v[28.81,30.15] s[45.89,46.94]
3.20s f[  0.31,  0.34] a[8.25,9.55] v[29.63,31.11] s[48.81,50.00]
3.30s f[  0.33,  0.36] a[8.16,9.53] v[30.45,32.06] s[51.81,53.16]
3.40s f[  0.35,  0.39] a[8.07,9.51] v[31.26,33.01] s[54.90,56.41]

Jak widać powyżej, po spadnieciu o 2-3 m (0.7-0.8 sek) różnica jest nierejestrowalna, dopiero po opadnieciu o ok. 50 m (3,3 sek) zarejestujemy różnicę gołym okiem.

Natomiast gdyby cieższa butelka (złośliwie ;-), obróciła się bokiem do kierunku lotu (przekrój wzrasta do ok. 0.025 m2) a lżejsza leciała w pozycji opływowej to wynik jest zależny od współczynika oporu i z tą dokładnościa obliczeń (szacowanie współczynników oporu) nie można go prawidłowo przewidzieć.