Computing work done by a force along a path

23. Computing work done by a force along a path#

Now that we know how to perform numerical differentiation and integration, we can begin applying the techniques to solve engineering problems. We start with a classic application: computing work.

Solar pressure#

Our sun continuously radiates energy in the form of electromagnetic waves. When these waves are reflected off of or absorbed by a surface, pressure is applied. These infinitesimal forces cause a small acceleration, which, in the case of a satellite in orbit, modifies its path.

In this example, we’ll compute the work done by the force due to solar pressure acting on a satellite traveling in an elliptical orbit around the sun.

The force due to solar radiation can be written as

\[ \vec{F} = -\dfrac{1}{c} GA \vec{s} \]

where:

  • \(c\) is the speed of light, \(3 \times 10^8\) m s\(^{-1}\).

  • \(G\) is the incident solar power flux, \(G = 1350\) Watt m\(^{-2}\).

  • \(A\) is the area of the spacecraft normal to the sun. The area is assumed to be constant \(A = 8\) m\(^2\).

  • \(\vec{s}\) is the unit vector directed from the spacecraft toward the sun.

The orbit is planar and is given parametrically by the following equation:

\[ r(t) = 1.5 \times 10^{11} \cos(t)\ \vec{i} + 1.3 \times 10^{11} \sin(t)\ \vec{j} \quad \text{with}\ 0 \le t \le 2\pi \]

Find the work done by the force of solar radiation on the spacecraft over the following interval: \(0 \le t \le \dfrac{\pi}{2} \). Use the time step of \(dt = 0.001\).

Procedure#

To calcuate the work done, we apply the following formula:

\[ W = \int \vec{F} \cdot \mathrm{d}\vec{r} = \int \vec{F} \cdot \vec{v}\ \mathrm{d}t = \int |\vec{F}| |\vec{v}| \cos\theta\ \mathrm{d}t, \]

where \(\theta\) is the angle between the velocity vector and the force. Each component of the integrand is computed separately. Then the integral is computed using the trapezoid() command.

import numpy as np
from scipy.integrate import trapezoid

# constants
G = 1350
A = 8
c = 3e8
dt = 0.001

# orbit
t = np.arange(0, np.pi/2, dt)
a = 1.5e11
b = 1.3e11
x = a * np.cos(t)
y = b * np.sin(t)

# velocity
vx = np.gradient(x, dt)
vy = np.gradient(y, dt)
speed = np.sqrt(vx**2 + vy**2)

# acceleration
ax = np.gradient(vx, dt)
ay = np.gradient(vy, dt)
accel = np.sqrt(ax**2 + ay**2)

# force mag and angle between v and F (use a as proxy)
F = -1/c * G * A
costheta = (vx * ax + vy * ay) / speed / accel

# work done
W = trapezoid(np.abs(F) * speed * costheta, t)
print(f"The work done is {W:.2f} J")

The work done is 720001.14 J

Note

Due to slight differences in implementation, particularly with handling boundary values, the result from trapezoid() in Python will differ from the result obtained using trapz() in MATLAB. 🤷🏼‍♂️

Contents