As already stated I am not familiar with a simpler way than the standard derivation. I don’t even think it is possible to derive it in a easier manner: For Stokes’ formula it is necessary to find expressions for pressure and viscous friction that to my knowledge can only be deducted from the Stokes’ equation and summing up their contribution to their drag force by integration. I would try to explain them the context and leave the mathematical formulas to the interested students.
I have not much idea about their background but I would structure my explanation as follows:
I would try to make them understand that the Stokes’ formula basically comes from integrating the stresses along the contour of a sphere. Theoretically if you know the flow field this approach can be used to determine resulting forces on any structure. Finding the stress distribution for a flow field analytically poses the biggest challenge and so far has only been shown for very low Reynolds number flow, where friction dominates, and the most symmetric of all shapes in 3D, a simple sphere. I think the most interesting part is understand how this is done and not precisely every single step in detail. If they follow you until my explanation starts to introduce the Navier-Stokes equations that is already perfectly fine. In the end most of the rest is just math that likely exceeds their current state of knowledge.
Fluid flow: smooth or eddies? - Reynolds number
In a fluid there are two dominating effects, one is inertia and the other is friction. Inertia characterises the desire of a fluid to keep its current state, its current velocity and direction of flow whereas friction basically determines if nearby fluid elements feel the state of the other, how much influence they have on each other, basically the “cohesion” between the fluid which counters the inertia.
In fluid dynamics you try to characterise the behaviour of a fluid by means of dimensionless numbers that are indifferent to the precise magnitude of single parameters. If all relevant dimensionless parameters in a model are similar its behaviour is also similar. This allows engineers and researchers to experiment with small-scale models instead of having to build full-scale airplanes and ships.
The one that characterises if a fluid is dominated by inertia or friction is the Reynolds number.
$$ Re := \frac{ U L }{ \nu } = \frac{\rho \, U \, L}{\mu} = \overbrace{\rho A U^2}^\text{inertia} \cdot \underbrace{\frac{L}{\mu A U}}_\text{viscous forces} \propto \frac{\text{inertial forces}}{\text{viscous damping forces}}.$$
For high Reynolds numbers inertia dominates far from the walls but near the walls the viscous friction takes over. As a result the flow forms eddies near the walls that are carried downstream and result in chaotic motion, turbulence. On the other hand for low Reynolds numbers viscosity (friction) prevails across the entire flow field: The flow is smooth, flows in layers (termed laminar), and even time-reversible (video of laminar demixing)
As you can see there are three parameters that influence this behaviour: The macroscopic velocity, the length scale of the problem and the viscosity of the medium. If we move a spoon with the same velocity through water and oil the outcome will be very differently: Water has a lower viscosity and thus the flow is characterised by a higher Reynolds number while the flow in oil will be characterised by a lower Reynolds number.
Educated guess and dimensional analysis
The first thing we can attempt to do is trying to make an educated guess what the relationship between drag force and the different parameters might look like using dimensional analysis. The unit of force is Newton
$$ [F_D] = N = \frac{kg \, m}{s^2}.$$
A correct relationship must yield the same units. We might assume that parameters influencing the outcome will be the diameter of the sphere $[D] = m$, the velocity $[U] = \frac{m}{s}$, the viscosity $[\nu] = \frac{m^2}{s}$ and the density $[\rho] = \frac{kg}{m^3}$. Larger values of all of the variables states above should result in higher drag force. Thus it seems plausible that such a relationship might take the form
$$ F_D \propto \nu \rho U D = \mu U D.$$
Yet we can already see that this is only one potential correlation. For higher Reynolds numbers actually the correlation looks like
$$ F_D \propto \rho U^2 A $$
where $A$ is the projected area with units $[A] = m^2$. While we could try to conduct a set of controlled experiments to verify our claims, we will follow a more formal way of deriving the drag force that is valid for any flow regime of continuum flows, the integration of surface stresses.
Forces on structures: Integration of stresses
A fluid flowing around an object exerts forces on its surface, which change with the flow field depending on the precise position. On one hand the pressure acts perpendicular to the surface and on the other hand the flowing fluid exerts a tangential friction force. Also the resulting force theoretically may have two components, one parallel to the flow, the drag and one perpendicular to it, the lift. If the object is designed properly and the flow conditions are chosen accordingly (even a flat plate is able to fly if you keep it at a certain angle to the flow) such that the pressure distribution results in a different pressure on the two sides of the object the lift can be very large such that it may make an object fly like the wings of an airplane (low pressure above and high pressure below) or help with traction in race cars such as spoilers in cars (upside-down wings).
In order to evaluate the resulting force we have to sum up all the tiny forces per surface element, the stresses, by means of integration (I have created an animated version of this figure that you can access here).
Our goal is now to find an analytical description for the drag. As we can imagine for highly turbulent (high Reynolds number) flow this is almost impossible to obtain as the flow will be chaotic, with no true steady state and instead depend on forming eddies. On the other side, for very low Reynolds numbers where the flow is smooth and follows the shape without separation we might be more lucky. In this flow regime we can use the rotational symmetry of the flow ($\Phi$ in this case is an angle around the horizontal rotation axis)
$$
u_{\Phi}=0, \quad \frac{\partial}{\partial \Phi}=0
$$
and the corresponding forces it suffices to integrate the forces over tiny cylinders surfaces. This can be done by looking at a cylindrical segment whose area is given by a circle with radius $R*\sin(\theta)$ and a height equal to the differential arc length $R d\theta$ (here a YouTube video about it and see this post on Archimedes' Hat-Box Theorem).
$$dA = \underbrace{2 R \sin(\theta) \pi}_\text{circle} \overbrace{R d\theta}^\text{height}$$
We are only interested in the drag and thus only this component is evaluated. In the considered flow regime the flow is symmetric anyways and thus the components in direction of the lift force would compensate. Note that this is not necessarily the case for all flow regimes: There exists a certain flow regime where the flow behind a cylinder forms an oscillating instability, the van Karman vortex street which leads to an oscillating lift force.
The resulting contributions to the drag force are given by:
$$ F_D = - \int\limits_{0}^{\pi} \sigma_{r \theta} |_{r=R} \sin(\theta) dA + \int\limits_{0}^{\pi} \sigma_{r r} |_{r=R} \cos(\theta) dA $$
$$ F_D = - 2 \pi R^2 \int\limits_{0}^{\pi} \sigma_{r \theta} |_{r=R} \sin^2(\theta) d\theta + 2 \pi R^2 \int\limits_{0}^{\pi} \sigma_{r r} |_{r=R} \sin(\theta)\cos(\theta) d\theta $$
We would have to evaluate these integrals now. The question is: What does the distribution of pressure look like? Well, the distribution can take a very complex form - it depends on how energy is converted between pressure and velocity and how much of it is dissipated (converted into heat) by friction. This depends on the flow parameters as well as the precise geometry. We will need some more fundamental correlations. This is where things get tricky and highly mathematical.
Navier-Stokes-Fouriers equations: The PDEs governing fluid flow
Physics are generally described by partial differential equations or PDEs for short. PDEs are equations that not only balance a variable but also its changes. E.g. an equation that does not only depend on position but also velocity (which is the change in position) and acceleration (which is the change in velocity). For instance take a simple object and air resistance which depends on velocity.
$$m \vec a = \vec F_D$$
If this force $\vec F_D$ was a constant value this equation could be solved very easily by double integration. Generally the drag for though depends on the velocity $\vec F_D = - k |\vec u| \vec u$ though and thus acceleration $\vec a = \frac{d \vec u}{dt}$ and velocity $u$ are coupled which makes them very complicated to solve analytically. For this reason one often turns to numerical solutions.
In case of fluid mechanics the governing equations are the Navier-Stokes-Fourier equations
$$\frac{\partial \rho}{\partial t} + \sum\limits_{j \in \mathcal{D}} \frac{\partial (\rho u_j )}{\partial x_j }=0 \tag{1}\label{1}$$
$$\frac{\partial (\rho u_i )}{\partial t} + \sum\limits_{j \in \mathcal{D}}\frac{\partial (\rho u_i u_j )}{\partial x_j} = \sum\limits_{j \in \mathcal{D}} \frac{\partial \sigma_{ij}}{\partial x_j } + \rho g_i \tag{2}\label{2}$$
$$\frac{\partial (\rho e)}{\partial t} + \sum\limits_{j \in \mathcal{D}} \frac{\partial (\rho u_j e)}{\partial x_j} = - \sum\limits_{j \in \mathcal{D}} \frac{\partial q_j}{\partial x_j} + \sum\limits_{i, j \in \mathcal{D}} \frac{\partial (\sigma _{ij} u_i)}{\partial x_j} + \sum\limits_{j \in \mathcal{D}} \rho u_j g_j \tag{3}\label{3}$$
where the total energy is given by the combination of internal $e_{in}$ and macroscopic energy $e := e_{in} + \sum\limits_{j \in \mathcal{D}} \frac{u_j u_j}{2}$, the local heat flux $q_i$ is generally assumed to be proportional to the gradient of the transported quantity, in this case the temperature, according to Fourier's law
$$q_i = - k \frac{\partial T}{\partial x_i}$$
and the stress tensor $\sigma_{ij}$ is given for a Newtonian fluid (detailed derivation here) by
$$\sigma_{ij} = - p \delta_{ij} + 2 \mu S_{ij} - \frac{2}{3} \mu \sum\limits_{k \in \mathcal{D}} S_{kk} \delta_{ij}. \tag{4}\label{4}$$
They basically describe the conservation of mass \eqref{1}, momentum \eqref{2} and energy \eqref{3} for a continuum fluid. The transport through temporal changes and advection (left-hand side) must equal the sources and sinks on the right-hand side (further explanation here). As you can see almost every term describes a change of a certain quantity (derivatives) which makes the system highly coupled.
Incompressible Navier-Stokes-equations
A first simplification is the assumption of an incompressible fluid which basically implies that density is a constant throughout the flow field. We furthermore silently assume that body forces can be neglected and the temperature does not change significantly throughout the flow field. As a consequence we only need the continuity and momentum equation which degenerate to:
$$\sum\limits_{j \in \mathcal{D}} \frac{\partial u_j}{\partial x_j} = 0,$$
$$\frac{\partial u_i}{\partial t} + \sum\limits_{j \in \mathcal{D}} u_j \frac{\partial u_i}{\partial x_j} = - \frac{1}{\rho} \frac{\partial p}{\partial x_i } + \frac{1}{\rho} \sum\limits_{j \in \mathcal{D}} \frac{\partial \tau_{ij}}{\partial x_j }.$$
and furthermore the last term in the stress tensor vanishes in equation \eqref{4} (for more details see).
Stokes’ equations: The PDEs of creepying flows ($Re \ll 1$)
Now our goal is to reduce complexity by checking if we can neglect some terms in these equations as others are more important and will dominate the flow for very low Reynolds number anyways. For very low Reynolds number $Re \ll 1$ clearly the viscosity will dominate over the inertia. In order to determine precisely which terms can be neglected we convert the equations to their dimensionless form by dividing through characteristic measures as follows:
$$x_i^*=\frac{x_i}{L}, \phantom{abc} u_i^*=\frac{u_i}{U}, \phantom{abc} \rho^*=\frac{\rho}{\rho_0}, \phantom{abc} t^*=\frac{t}{\frac{L}{U}},\phantom{abc} p^*=\frac{p}{\rho_0 \nu \frac{U}{L}} \phantom{ab}$$
This results in the two dimensionless equations
$$\sum\limits_{j \in \mathcal{D}} \frac{\partial u_j^*}{\partial x_j^* }=0$$
$$\underbrace{\overbrace{Re}^{\ll 1} \left( \frac{\partial u_i^*}{\partial t^*} + \sum\limits_{j \in \mathcal{D}} u_j^* \frac{\partial u_i^*}{\partial x_j^*} \right)}_{\ll 1} = - \frac{\partial p^*}{ \partial x_i^* } + \sum\limits_{j \in \mathcal{D}} \frac{\partial \tau_{ij}^*}{\partial x_j^* }$$
where the last one can be rewritten with equation \eqref{4} and neglecting the terms on the left-hand side due to their magnitude to
$$\frac{\partial p}{\partial x_i} = \mu \sum\limits_{j \in \mathcal{D}} \frac{\partial^2 u_i}{\partial x_j^2}.$$
This set of equations is often termed Stokes' equations and can be written symbolically as
$$\vec \nabla \cdot \vec u = 0,$$
$$\vec \nabla p = \mu \vec \nabla^2 \vec u = \mu \vec \Delta \vec u.$$
Now it is possible to find a closed descriptions for the pressure and the velocity by solving this much simpler system of partial differential equations.
Vector identities and rotational symmetry
We can use the vector identity
$$\vec \nabla ^2 \vec u = \vec \nabla ( \underbrace{ \vec \nabla \cdot \vec u }_{= 0}) - \vec \nabla \times ( \vec \nabla \times \vec u)$$
to simplify the Stokes' momentum equation to
$$\vec \nabla p = - \mu \vec \nabla \times ( \vec \nabla \times \vec u).$$
Now we can take the cross-product of this equation and apply the vector identity
$$\vec \nabla \times \vec \nabla p = 0.$$
in order to eliminate the pressure and obtain the linear equation
$$\vec \nabla \times [ \vec \nabla \times ( \vec \nabla \times \vec u) ].$$
Further we will now switch to spherical coordinates as it is more convenient for a rotational symmetric problem
$$\vec u = \left(\begin{array}{c}
u_r \\
u_\Theta \\
u_\Phi
\end{array}\right)$$
The corresponding operators take the following form
$$\require{cancel}
\vec \nabla \cdot \vec u =\frac{1}{r^{2}} \frac{\partial}{\partial r}\left(r^{2} u_{r}\right)+\frac{1}{r \sin \Theta} \frac{\partial}{\partial \Theta}\left(u_{\Theta} \sin \Theta\right)+\cancel{\frac{1}{r \sin \Theta} \frac{\partial u_{\Phi}}{\partial \Phi}},$$
$$\require{cancel}
\vec \nabla \times \vec u =\left(\begin{array}{c}
\cancel{\frac{1}{r \sin \Theta} \frac{\partial}{\partial \Theta}\left(u_{\Phi} \sin \Theta\right)}-\cancel{\frac{1}{r \sin \Theta} \frac{\partial u_{\Theta}}{\partial \Phi}} \\
\cancel{\frac{1}{r \sin \Theta} \frac{\partial u_{r}}{\partial \Phi}}-\cancel{\frac{1}{r} \frac{\partial}{\partial r}\left(r u_{\Phi}\right)} \\
\frac{1}{r} \frac{\partial}{\partial r}\left(r u_{\Theta}\right)-\frac{1}{r} \frac{\partial u_{r}}{\partial \Theta}
\end{array}\right).$$
Stokes' stream function and solving the PDE
The derivation from here on can be found in a similar way in this NYU lecture notes, section 7.3. We can now introduce Stokes' stream function for rotational symmetric bodies
$$\frac{\partial \Psi}{\partial \Theta}:=u_{r} r^{2} \sin \Theta, \quad \frac{\partial \Psi}{\partial r}:=-u_{\Theta} r \sin \Theta.$$
This allows us to rewrite the Stokes' momentum equation as follows:
$$\vec \nabla \times \vec u =\left(\begin{array}{c}
0 \\
0 \\
\frac{1}{r} \frac{\partial}{\partial r}\left(r u_{\Theta}\right)-\frac{1}{r} \frac{\partial u_{r}}{\partial \Theta}
\end{array}\right)=
\left(\begin{array}{c}
0 \\
0 \\
-\frac{1}{r sin \Theta} \underbrace{ \left[ \frac{\partial^2 \Psi}{\partial r^2} + \frac{\sin \Theta}{r^2} \frac{\partial}{\partial\Theta} \left( \frac{1}{\sin \Theta} \frac{\partial \Psi}{\partial \Theta} \right) \right] }_{:= \mathcal{L}(\Psi)}
\end{array}\right),$$
$$\vec \nabla \times(\vec \nabla \times \vec u)=\left(\begin{array}{c}
- \frac{1}{r^2 \sin \Theta} \frac{\partial [\mathcal{L}(\Psi)]}{\partial \Theta} \\
\frac{1}{r \sin \Theta} \frac{\partial [ \mathcal{L}(\Psi)]}{\partial r} \\
0
\end{array}\right),$$
$$\vec \nabla \times[\vec \nabla \times(\vec \nabla \times \vec u)]=\left(\begin{array}{c}
0 \\
0 \\
\frac{1}{r \sin \theta} \underbrace{ \left\{\frac{\partial^{2}(\mathcal{L} \Psi)}{\partial r^{2}}+\frac{\sin \theta}{r^{2}} \frac{\partial}{\partial \Theta}\left[\frac{1}{\sin \theta} \frac{\partial(\mathcal{L} \Psi)}{\partial \Theta}\right]\right\} }_{= \mathcal{L}[\mathcal{L}(\Psi)]} \\
\end{array}\right).$$
This can also be written as
$$\vec \nabla \times[\vec \nabla \times(\vec \nabla \times \vec u)] =\frac{1}{r \sin \Theta} \mathcal{L}[\mathcal{L} (\Psi)] \vec e_{\Phi}
=\frac{1}{r \sin \Theta} \mathcal{L}^{2} \Psi \vec e_{\Phi}$$
where the operator $\mathcal{L}$ is given as
$$\mathcal{L}=\frac{\partial^{2}}{\partial r^{2}}+\frac{\sin \theta}{r^{2}} \frac{\partial}{\partial \Theta}\left(\frac{1}{\sin \Theta} \frac{\partial}{\partial \Theta}\right).$$
Additionally the flow has to fulfill the boundary conditions. The fluid velocity at the wall must be zero (no-slip condition)
$$u_r = 0: \frac{\partial \Psi}{\partial \Theta} = 0, \quad u_\Theta = 0: \frac{\partial \Psi}{\partial r} = 0$$
and in the far-field the velocity is given by the unperturbed velocity $U_\infty$
$$u_{r} =\frac{1}{r^{2} \sin \Theta} \frac{\partial \Psi}{\partial \Theta}=U_\infty \cos \Theta,$$
$$u_{\Theta} =-\frac{1}{r \sin \Theta} \frac{\partial \Psi}{\partial r}=-U_\infty \sin \Theta.$$
Solving each of the equations for $\Psi$ by means of integration we find
$$\Psi=U_\infty \frac{r^{2}}{2} \sin ^{2} \Theta+C \quad \forall \Theta, r \rightarrow \infty$$
for the solution indefinitely far from the cylinder. Thus, the resulting partial differential equation could be solved by a product ansatz of the form
$$\Psi(r, \Theta)=f(r) \sin ^{2} \Theta.$$
We insert this ansatz into the partial differential equation
$$\mathcal{L} (\Psi)=\underbrace{\left(f^{\prime \prime}-\frac{2}{r^{2}} f\right)}_{:= F(r)} \sin ^{2} \Theta,$$
$$\mathcal{L}^{2} (\Psi)=\mathcal{L}[\mathcal{L} (\Psi)]=\left(F^{\prime \prime}-\frac{2}{r^{2}} F\right) \sin ^{2} \Theta,$$
$$\mathcal{L}^{2} \Psi=0 \Leftrightarrow(\underbrace{F^{\prime \prime}-\frac{2}{r^{2}} F}_{= 0}) \underbrace{\sin ^{2} \Theta}_{=\neq 0}=0.$$
As the last term can't be assumed to vanish the differential equation that must be fulfilled is the simple Eulerian differential equation
$$F^{\prime \prime}-\frac{2}{r^{2}} F=0$$
that can be solved with the ansatz $F = C r^\lambda$ resulting in the algebraic equation
$$\cancel{r^2} \lambda (\lambda -1) \cancel{Cr^{\lambda-2}} - 2 \cancel{Cr^\lambda} = 0,$$
$$\lambda (\lambda - 1) - 2 = (\lambda - 2) (\lambda + 1) = 0,$$
which can be solved in a similar manner additionally considering the particular solution to
$$F = f^{\prime \prime}-\frac{2}{r^{2}} = A r^2 + \frac{B}{r}$$
which again results in
$$f(r)=\frac{A}{10} r^{4}-\frac{B}{2} r+C r^{2}+\frac{D}{r}.$$
Thus we find for the stream function $\Psi$ and the radial $u_r$ and tangential $u_\Theta$ velocities
$$\Psi=\frac{1}{4} U_\infty R^{2}\left(\frac{R}{r}-3 \frac{r}{R}+2 \frac{r^{2}}{R^{2}}\right) \sin ^{2} \Theta,$$
$$u_{r}=U_\infty \left(1+\frac{1}{2} \frac{R^{3}}{r^{3}}-\frac{3}{2} \frac{R}{r}\right) \cos \Theta,$$
$$u_{\Theta}=U_\infty \left(-1+\frac{1}{4} \frac{R^{3}}{r^{3}}+\frac{3}{4} \frac{R}{r}\right) \sin \Theta.$$
Finally we can determine the pressure by integrating the radial Stokes' momentum equation
$$\frac{\partial p}{\partial r} =\mu \frac{1}{r^{2} \sin \Theta} \frac{\partial[\mathcal{L} (\Psi)]}{\partial \Theta},$$
to
$$p =-\frac{3}{2} \frac{\mu U_{\infty} R}{r^{2}} \cos \theta+D$$
Considering the pressure in the far field $r \to \infty$ we finally yield
$$p = p_{\infty}-\frac{3}{2} \frac{\mu U_{\infty} R}{r^{2}} \cos \Theta.$$
The pressure in front is clearly the highest similar to the impulse you feel when sticking the hand out of a driving vehicle. Behind the sphere it is anti-symmetrically lower.
Integrating the force
With these distributions for pressure and velocity finally the stresses and therefore the drag force can be evaluated by integration
$$\left. \sigma_{rr} \right|_{r=R} = - \left. p \right|_{r=R} + \left. \underbrace{ 2 \mu \frac{\partial u_r}{\partial r} }_{\tau_{rr}} \right|_{r=R} = - p_{\infty} + \frac{3}{2} \frac{\mu U_{\infty}}{R} \cos \Theta,$$
$$\left. \sigma_{r \Theta} \right|_{r=R} = \left. \underbrace{ \mu \left( \frac{1}{r} \frac{\partial u_r}{\partial \Theta} + \frac{\partial u_\Theta}{\partial r} - \frac{u_\Theta}{r} \right)}_{\tau_{r \Theta}} \right|_{r=R} = - \frac{3}{2} \frac{\mu U_\infty}{R} \sin \Theta,$$
$$F_D = 3 \pi \mu U_{\infty} R \int\limits_{\Theta = 0}^{\pi} \underbrace{(\sin^2 \Theta + \cos^2 \Theta)}_{=1} \sin \Theta d\Theta = 3 \pi \mu U_{\infty} R \left. cos \overline \Theta \right|_{\overline \Theta = 0}^{\pi} = 6 \pi \mu U R = 3 \pi \mu U D$$
As can be seen this indeed takes the form that we predicted by smart dimensional analysis. Evaluating each contribution independently we find
$$F_D = \underbrace{ 4 \pi \mu U R}_\text{viscous contribution} + \underbrace{2 \pi \mu U R}_\text{pressure contribution}.$$
