Modelling assumptions in mechanics

Modelling in Mechanics

In virtually all mechanics problems some modelling assumptions need to be made in order to simplify the problems to a point where they can be analysed. Here we will consider two of the more common models that are used in the context of the Model Solutions problem - constant gravitational acceleration and neglection of air resistance - and look at how they affect the validity of the analysis.

Some of the mathematics here is quite involved (knowledge of vector notation and how to solve linear second order differential equations will be required) and also often quite long winded - this has intentionally been left the case however to demonstrate the importance modelling assumptions play in simplifying the mathematical analysis of even apparently simple physical situations.

Constant/variable gravitational acceleration

The universal gravitation law is:

$$\mathbf{F_g} = \frac{-GMm}{r^2} \mathbf{\hat{r}}$$

Where $\mathbf{F_g}$ is the force acting on a body of mass $m$ due to gravitational attraction to a second body of mass $M$ with vector from the centre of mass of the second body to that of the first $\mathbf{r} = r \mathbf{\hat{r}}$ (with $\mathbf{\hat{r}}$ being a unit vector) . $G$ is the gravitational constant and in SI units has value $6.67300 \times 10^{-11} \textrm{ m}{^3} \textrm{ kg}^{-1} \textrm{ s}^{-2}$.

The gravitational acceleration experienced by the body is therefore:

$$\mathbf{a_g} = \frac{-GM}{r^2} \mathbf{\hat{r}}$$

It can be seen that this acceleration is a function of the separation between the bodies, specifically there is an inverse square law relationship. However in most problems (based on Earth!) we assume a constant acceleration due to gravity, $g = 9.81 \textrm{ ms}^{-2}$. How can this be valid?

The assumption comes from the composition of the separation $r$ term for bodies in the Earth's gravitational field. For an object at an altitude $h$ above the Earth's surface, $r = h + R_E$ where $R_E$ is the radius of the Earth (approximately 6738 km). In general however $h \ll R_E$ and therefore it is valid to assume $r \approx R_E$ to a good degree of accuracy (for instance even for $h = 10 \textrm{ km}$, the cruise altitude of a jet airliner, the error in the acceleration associated with assuming $r = R_E$ is less than 0.2%).

Using the figures for $G$ and $R_E$ above and mass of Earth of $M_E = 5.9742 \times 10^{24} \textrm{ km}$ gives:

$$g = \frac{GM_E}{R_E^2} = 9.80 \textrm{ ms}^{-2}$$

(Note: The difference from the conventional value of $g = 9.81 \textrm{ ms}^{-2}$ arises from the variation in the Earth's radius with latitude due to it's non-spherical shape, with the standardised value using a slightly different $R_E$ from that quoted).

To show how assuming constant $g$ affects the analysis of the motion of objects, we will consider the motion of the shot put from the Model Solutions problem with and without the assumption. For now we will make the further modelling assumptions of a point mass with no air resistance to keep things (relatively!) simple.

Let the shot be projected at initial velocity $\textbf{v}(0) = p \textbf{i} + q \textbf{j}$

Case 1: Constant g

By symmetry arguments, the shot will take the same time to travel from its launch point to the apex of its trajectory, as from the apex to when it is again level with the point it was launched from. At the apex it will have instantaneous vertical speed $v = 0$ with the time it takes to reach there $\frac{1}{2}T$. $u$, $v$ and $a$ are known and $t$ is wanted therefore use constant acceleration formula $v = u + at$.

$$ \Rightarrow 0 = q + (-g) \frac{1}{2}T$$

Time taken to return to the same height as the launch point is therefore $T = \frac{2q}{g}$.

There is no net horizontal force on the shot throughout its flight and so its horizontal velocity remains constant at $p$. The horizontal range of the shot, $d$, is therefore:

$$d = pT = \frac{2pq}{g}$$

Case 2: Variable g

As established above:

$$a_g = \frac{\mathrm{d}^2 r}{\mathrm{d}t^2} = \frac{-GM_E}{r^2}$$

This is a non-linear second order ordinary differential equation and with the boundary conditions applicable in this case ($r(0) = R_E$ and $\dot{r}(0)=q$) it has no analytical solution in terms of elementary functions. To proceed we must therefore use a numerical technique to obtain an approximated solution.

Here we will use what is called Verlet integration, a numerical method often used in the computations involved in modelling motion.

The displacement $r$ of the particle is some function of time $t$. Considering $r$ at some time $t + \delta t$, we can form a Taylor series expansion:

$$r \left( t + \delta t \right) = r(t) + \frac{\mathrm{d}r(t)}{\mathrm{d}t} \delta t + \frac{\mathrm{d}^2 r(t)}{\mathrm{d}t^2} \frac{\delta t ^2}{2!} + \frac{\mathrm{d}^3 r(t)}{\mathrm{d}t^3} \frac{\delta t ^3}{3!} + O(\delta t ^4)$$

And similarly:

$$r \left( t - \delta t \right) = r(t) - \frac{\mathrm{d}r(t)}{\mathrm{d}t} \delta t + \frac{\mathrm{d}^2 r(t)}{\mathrm{d}t^2} \frac{\delta t ^2}{2!} - \frac{\mathrm{d}^3 r(t)}{\mathrm{d}t^3} \frac{\delta t ^3}{3!} + O(\delta t ^4)$$

Adding these two together gives:

$$r \left( t + \delta t \right) + r \left( t - \delta t \right) = 2r(t) + \frac{\mathrm{d}^2 r(t)}{\mathrm{d}t^2} \delta t ^2 + O(\delta t ^4)$$

Rearranging, substituting in the expression derived for acceleration and assuming $\delta t$ is small so terms in $\delta t^4$ and above can be ignored:

$$r \left( t + \delta t \right) = 2 r(t) - r(t - \delta t) - \frac{GM \delta t^2}{r(t)^2}$$

We now have an equation we can apply iteratively to compute how $r$ varies with $t$ given two starting values, so to find the range we continue iterating until the object returns to it's starting height. Although it would of course be possible to do this by hand, this would be extremely time consuming and prone to human error, with thousands of iterations likely to be needed. Such a task is however ideally suited to being carried out on a computer and this is how the figures that will be used for comparison below have been generated. Although not directly relevant to the article for those who are interested at how this might be done below is an example in Python.

    #Constants
    G = 6.673e-11
    M = 5.9742e+24
    R = 6.378e+06
    dt = 0.001
    #Function to calculate range of particle given initial
    #velocity v = pi + qj under constant gravity
    def calcRangeConst(p, q):
        g = (G*M)/(R*R)
        return 2*p*q/g
    #Function to calculate range of particle given initial
    #velocity v = pi + qj under variable gravity
    def calcRangeVar(p, q):
        t=dt
        hp = R
        h = hp + q*dt + 0.5*(-(G*M)/(R*R))*dt*dt
        while h>     =R:
            temp = h
            a = -(G*M)/(h*h)
            h = 2*h - hp + a*dt*dt
            hp = temp
            t+=dt
        return t*p

Comparing the results

Having now established methods for calculating ranges of a projected particle under both variable and constant gravitational acceleration, we can now identify the effect that assuming constant gravitational acceleration has. Below are calculated values for the range of a particle projected at various initial speeds under both constant and variable gravitational acceleration.

As suggested by the negligible change in gravitational acceleration over even very large height differences, the effect of modelling gravity as constant has virtually no effect on range. The diecrepancy increases with increasing vertical projection speed (with horizontal projection speed having no effect as would be expected), but even at quite significant speeds the percentage difference is tiny. It therefore seems fair to say that it is valid to model projectile motion as being under constant gravity in virtually all cases. The error introduced into results is virtually certain to be insignificant in comparison to those introduced by other necessary modelling assumptions. An equally important point is that from the working above it is evident that this model simplifies the range calculations massively, without this assumption it not in fact possible to solve the problem analytically at all.

Including air-resistance

The neglection of air resistance is a fairly major modelling assumption to make as drag effects can have a significant effect on an object's trajectory. However to accurately account for the effects of air resistance is very difficult as it has a complex dependence on the velocity of the object and of the properties of the fluid it is moving through, with only numerical rather than analytical solutions possible.

For low velocity motion through a fluid such as air, a very basic model of air resistance where the magnitude of the resistive force is proportional to the relative velocity between the object and fluid and is always directed to oppose motion, can be used to give a general indication of it's effects. This is called Stoke's drag, and the relevant equation is:

$$\mathbf{F_{drag}} = - b \mathbf{v}\,,$$

where $\mathbf{F_{drag}}$ is the drag force experienced by an object moving with velocity $\mathbf{v}$ relative to the fluid it is moving through. The proportionality constant $b$ depends on the viscosity of the fluid and geometry of the object.

In general air resistance will decrease both the horizontal range and time in air of the shot put, with the level of the effect depending on the value of the proportionality constant $b$. Below the equations of motion for a shot-put thrown both with and without air resistance considered are derived to demonstrate this.

As before, let the shot be projected at initial velocity $\textbf{v}(0) = p \textbf{i} + q \textbf{j}$ from position $\mathbf{r}(0) = 0 \mathbf{i} + 0 \mathbf{j}$

Case 1: No air resistance

Applying Newton's Second Law to the shot in the vertical direction:

$$m\ddot{y} = -mg \Rightarrow \ddot{y} = -g$$

Integrating with respect to time twice:

$$\dot{y} = -gt + A$$ $$\ddot{y} = -\frac{1}{2}gt^2 + At + B$$

Applying initial conditions to find constants:

$$y(0) = 0 \Rightarrow 0 = -g(0) + A \Rightarrow A = 0$$ $$\dot{y}(q) = \frac{1}{2}g(0) + A(0) + B \Rightarrow B = q$$ $$\therefore y = qt - \frac{1}{2}gt^2$$

There are no forces acting on the shot in the horizontal direction, and so there is no horizontal acceleration. The horizontal speed therefore remains constant.

$$\therefore x = pt$$

The overall equation of motion is therefore:

$$\textbf{r}(t) = (pt) \textbf{i} + (qt - \frac{1}{2}gt^2) \textbf{j}$$

To calculate the range of the shot-put, the time $T$ at which the shot returns to ground level needs to be calculated

$$\Rightarrow qT -\frac{1}{2}gT^2 = 0 \Rightarrow T = \frac{2q}{g}$$

Substituting this into the equation for the horizontal motion, give the range:

$$D = pT = \frac{2pq}{g}$$

Case 2: Including air resistance

Applying Newton's Second Law to the shot put in the vertical direction:

$$m\ddot{y} = -mg - b\dot{y} \Rightarrow \ddot{y} + \frac{b}{m}\dot{y} = -g$$

This is a inhomogeneous linear second order differential equation in $y$ and can be solved using the complementary function / particular integral technique.

Complementary function:

$$\lambda^2 + \frac{b}{m}\lambda = 0 \Rightarrow \lambda\left(\lambda + \frac{b}{m}\right) = 0 \Rightarrow \lambda = 0 \textrm{ or } -\frac{b}{m}$$ $$\Rightarrow y_{CF} = Ae^{0t} + Be^{-\frac{b}{m}t} = A + Be^{-\frac{b}{m}t}$$

Particular integral:

$$y_{PI} = Ct \Rightarrow \dot{y}_{PI} = C ;\ \ddot{y}_{PI} = 0$$ $$\Rightarrow 0 + \frac{b}{m}C = -g \Rightarrow C = - \frac{mg}{b}$$

Superposition $y = y_{CF} + y_{PI}$:

$$y = A + Be^{-\frac{b}{m}t} - \frac{mg}{b}t$$ $$\dot{y} = \frac{-bB}{m}e^{-\frac{b}{m}t} - \frac{mg}{b}$$

Applying initial conditions to find constants:

$$y(0) = 0 \Rightarrow 0 = A + B(1) - \frac{mg}{b}(0) \Rightarrow A = -B$$ $$\dot{y}(0) = q \Rightarrow q = \frac{-bB}{m}(1) - \frac{mg}{b}$$ $$B = \frac{-m}{b} \left( \frac{mg}{b} + q \right) \Rightarrow A = \frac{m}{b} \left( \frac{mg}{b} + q \right)$$ $$\therefore y = \frac{m}{b} \left( \left( \frac{mg}{b} + q \right) \left( 1 - e^{-\frac{b}{m}t} \right) - gt \right)$$

Applying Newton's Second Law to the shot-put in the horizontal direction:

$$m\ddot{x} = -b\dot{x} \Rightarrow \ddot{x} + \frac{b}{m}\dot{x} = 0$$

This is again a linear second-order equation, this time homogeneous with a solution the same as that to the complementary function previously:

$$x = D + Ee^{-\frac{b}{m}t}$$ $$\dot{x} = -\frac{bE}{m}e^{-\frac{b}{m}t}$$

Applying initial conditions to find constants:

$$x(0) = 0 \Rightarrow 0 = D + E(1) \Rightarrow D = -E$$ $$\dot{x}(0) = p \Rightarrow p = \frac{-bE}{m}(1)$$ $$E = -\frac{mp}{b} \Rightarrow D = \frac{mp}{b}$$ $$\therefore x = \frac{mp}{b} \left( 1 - e^{-\frac{b}{m}t} \right)$$

The overall equation of motion is therefore:

$$\textbf{r}(t) = \frac{mp}{b} \left( 1 - e^{-\frac{b}{m}t} \right) \textbf{i} + \frac{m}{b} \left( \left( \frac{mg}{b} + q \right) \left( 1 - e^{-\frac{b}{m}t} \right) - gt \right) \textbf{j}$$

To find the time in air here, the time T when the vertical height is zero again needs to be found:

$$\Rightarrow \frac{m}{b} \left( \left( \frac{mg}{b} + q \right) \left( 1 - e^{-\frac{b}{m}T} \right) - gT \right) = 0$$ $$\Rightarrow \left( \frac{mg}{b} + q \right) \left( 1 - e^{-\frac{b}{m}T} \right) = gT$$ $$\Rightarrow 1 - e^{-\frac{b}{m}T} = \frac{gb}{mg + qb}T$$

This equation can not be solved analytically. However as this model is only valid for low drags, we can assume $\frac{b}{m}T$ is small. Therefore using the power series expansion for the exponential function, excluding terms of order four and higher:

$$\begin{align*}1 - \left( 1 + \left( \frac{-b}{m}T \right) + \frac{1}{2!}\left( \frac{-b}{m}T \right)^2 + \frac{1}{3!}\left( \frac{-b}{m}T \right)^3 \right) &= \frac{gb}{mg + qb}T\\ \frac{b}{m}T -\frac{b^2}{2m^2}T^2 + \frac{b^3}{6m^3}T^3 &= \frac{gb}{mg + qb}T\\ \frac{b^3}{6m^3}T^3 -\frac{b^2}{2m^2}T^2 + \left(\frac{b}{m} - \frac{gb}{mg + qb} \right) T &= 0\\ \frac{b}{m}T \left( \frac{b^2}{6m^2}T^2 - \frac{b}{2m}T + 1 - \frac{mg}{mg + qb} \right) &= 0 \\ \frac{b}{m}T \left( \frac{b^2}{6m^2}T^2 - \frac{b}{2m}T + \frac{mg + qb}{mg + qb} - \frac{mg}{mg + qb} \right) &= 0\\ \frac{b}{m}T \left( \frac{b^2}{6m^2}T^2 - \frac{b}{2m}T + \frac{qb}{mg + qb} \right) &= 0\\ \therefore T = 0 \textrm { (as expected) or } \frac{b^2}{6m^2}T^2 - \frac{b}{2m}T + \frac{qb}{mg + qb} &= 0\end{align*}$$$$\begin{align*}T &= \frac{3m^2}{b^2} \left(\frac{b}{2m} \pm \sqrt{\frac{b^2}{4m^2} - 4 \frac{b^2}{6m^2} \frac{qb}{mg + qb}}\right)\\T &= \frac{3m}{2b} \left(1 \pm \sqrt{1 - \frac{8qb}{3(mg + qb)}} \right)\\ T &= \frac{3m}{2b} \left(1 \pm \sqrt{\frac{3mg - 5qb}{3mg + 3qb}} \right)\end{align*}$$

Only the lower of these roots is a valid approximation (as $T$ increases the series approximation diverges from the actual result), i.e.

$$T = \frac{3m}{2b} \left(1 - \sqrt{\frac{3mg - 5qb}{3mg + 3qb}} \right)$$

To find the range $D$ of the projectile this time value should be substituted into the equation for the horizontal coordinate, this giving the horizontal distance covered at the point it returns to the ground:

$$D = \frac{mp}{b} \left( 1 - e^{-\frac{b}{m}T} \right)$$

Comparing the results

What is immediately obvious is that the analysis without air resistance is much shorter than when it is included!

Tabulated below are the calculated ranges of a projectile with various initial speeds using the equations derived above. Values are given for the range with no drag effect included, a low drag effect ($\frac{b}{m}=0.05$) and a high drag effect ($\frac{b}{m}=0.2$). The percentages in brackets indicate the decrease from the range with no drag effect at that projection velocity.

As we would intuitively expect, as the projectile's initial speed increases so does the discrepancy between the results with and without the drag effect included. The decreased range at higher $\frac{b}{m}$ values also fits with everyday experience: the more viscous the fluid being travelled through (corresponding to a higher $b$ value) and the lower the mass of the projectile, the greater the effect of drag - for instance a golf ball travelling through air (low $\frac{b}{m}$) in comparison to a table tennis ball travelling through water (high $\frac{b}{m}$).

Though the model used for drag here was a gross oversimplification of reality, it still gives sufficient evidence to show neglection of drag effects / air resistances in modelling projectile motion has a very significant affect on the results obtained even at relatively low velocities. Conversely however, it can be seen that even with such a simplistic model for air resistance, the maths involved in modelling projectile motion quickly becomes very difficult and time-consuming. The methods required to more accurately model drag effects are even more complex by several orders of magnitude and even these give only an approximation for idealised situations.

This conflict between accuracy and complexity is evident in nearly all modelling, the inverse relationship between the two meaning a compromise has to made between the desired accuracy of the results and time and effort available to obtain them. It is therefore always important to consider the effect of any modelling assumptions in the context of the intended end-use of the analysis, and to keep in mind any potential discrepancies in results when using them.