Hover over the menu bar to pick a physics animation.

Saved Force Data Vs Sphere Speed:
178.86,197.62,218.00,241.15,263.29
Saved Cd Data Vs Sphere Speed:
2.87,2.55,2.38,2.25,2.15

Drag on a Large Sphere Due to a Digital Gas

I have decided to change the way that I do this computation. I will apply a constant downward speed to the sphere and compute the gas forces needed to maintain that speed. This will be much better controlled than using a gravity force to induce the speed. Also, the sphere has to be a "rough" sphere or we would never be able to develop the laminar flow streamlines around it. For a slow moving sphere and with significant attractive forces between the gas particles we expect a relation called Stokes Law for the force, `F` on a sphere moving through the fluid .

`F=6pi a eta v \ \ \ (1)`

where `a` is the radius of the sphere, `eta` is the viscosity coefficient, and `v` is the speed of the sphere. In this case the fluid will be a digital gas that has no attractive forces between its particles. Using Stokes Law, from the speed and the applied resulting force of the gas on the sphere we can compute the viscosity coefficient of the digital gas.

A separate way of computing the force that a digital gas exerts on a fast moving sphere is to use the standard `v^2` dependent expression for drag force:

`F=1/2 C_d rho v^2 A \ \ \ (2)`

where `C_d` is the drag coefficient, `rho` is the gas density, `v` is the speed of the sphere and A is the frontal area of the sphere, `A=pi r_S^2` where `r_S` is the sphere radius. The literature, see Force for High Reynolds Number, has it that the sphere shape has an approximate drag coefficient `C_d=0.5`. See also Drag of a Golf Ball. If we know the gas drag force on the sphere , we can use equation 2 to compute the `C_d` for our digital gas.

`C_d= 2F/( rho v^2 A) \ \ \ (2a)`

where `v` is the constant velocity of the sphere resulting in the gas drag force `F`. As the references note, `C_d` is expected to change with sphere velocity. It's important to understand that the reference values of `C_d` refer to a sphere whose surface is rough compared to the size of the particles hitting it. Without the roughness the laminar flow induced streamlines around the sphere would never develop. For this animation the roughness will be a corrugation of radius variations that depend on the latitude angle `theta` on the sphere.

Computing the Corrugated Sphere Boundary and its Surface Normals

The boundary is quite simple to compute in terms of the sphere mean radius, `r_S`, the amplitude, `A` of the sinusoidal corrugations, and the number of corrugation periods, `n_P`. If we let `theta` be the latitude coordinate and `phi` be the longitude coordinate, then the radius of a point on the sphere will be

`r_(pt)=r_S+A cos(n_P theta) \ \ \ (3)`

where `theta` varies from `-pi/2` to `pi/2` and is independent of the value of `phi`. The (x,y,z) coordinates of this point will be:

`x=r_(pt)cos theta cos phi \ \ \ (4a)`
`y=r_(pt)cos theta sin phi \ \ \ (4b) `
`z=r_(pt)sin theta \ \ \ (4c)`

If we know (x,y,z) instead of `theta` and `phi` we can obtain `theta` and `phi` from the following equations.

`sin phi= y/sqrt(x^2+y^2) \ \ \ (5a)`
`cos phi= x/sqrt(x^2+y^2) \ \ \ (5b)`

`sin theta = z/sqrt(x^2+y^2+z^2) \ \ \ (6a)`
`cos theta = (x/cos phi)/sqrt(x^2+y^2+z^2) \ \ \ (6b)`

For a three dimensional corrugated sphere `theta` is limited to the range `-pi/2 to pi/2` so `cos theta` isalways greater than 0.

The surface normals depend on the derivative of `r_(pt)` with respect to `theta`.

`slope=tan psi=(dr_(pt))/(d(r_S theta))=-n_P A/r_S sin(n_P theta) \ \ \ (7)`

This expression is the slope of the same corrugation centered on the x axis. The sine and cosine of the slope angle are expressed simply as:

`cos psi=1/(sqrt(1+tan^2 psi) \ \ \ (8a)`
`sin psi=tan psi/(sqrt(1+tan^2 psi)\ \ \ (8b)`

where again `cos psi` is always greater than zero.

We need to transform this slope so it relates to the three unit normals of the corrugagted sphere at any point (x,y,z). Without explicit proof, I now state the components of the unit vector that is normal to the corrugated sphere in the form:

`bb hat "n"_r=n_x bb hat "x"+n_y bb hat "y"+n_z bb hat "z" \ \ \ (9)`

In terms of `theta, phi, and psi` this can be written:

`bb hat "n"_r=cos(theta-psi) cos phi\ bb hat "x"+cos(theta-psi)sin phi \ bb hat "y"+ sin(theta-psi) \ bb hat "z" \ \ \ (10)`

where using trigonometric identities we can write

`cos (theta-psi)=cos theta cos psi+sin theta sin psi \ \ \ (11a)`
`sin (theta-psi)=sin theta cos psi-cos theta sin psi) \ \ \ (11b)`

Note that expressions for `cos theta, sin theta, cos psi and sin psi` in terms of `(x,y,z)` are provided in the equations above.

Equation 10 is to be contrasted with the similar equation for a sphere without corrugation:

`bb hat "n"_r=costheta cos phi\ bb hat "x"+cos thetasin phi \ bb hat "y"+ sin theta \ bb hat "z" \ \ \ (12)`

Note that the normal vector in both equations (10) and (12) have unit magnitude which can be confirmed by taking the dot (inner) product of each vector with itself e.g.

`bb hat "n"_r cdot bb hat "n"_r=cos^2theta cos^2phi+cos^2theta sin^2phi + sin^2theta =1\ \ \ (13)`

Computing the Gas Drag Force on the Corrugated Sphere Boundary

Having obtained the surface boundary radius `r` and the normal unit vector `bb hat "n"_r`, the expression for the new particle velocity, `bb "v"_a`, vector after a collision of a particle with the sphere is

`bb "v"_(a)=bb "v"_(b) -2 (bb hat "n"_r cdot [bb "v"_b-bb "v"_S)] bb hat "n"_r \ \ \ (14)`

where `bb "v"_b` represents the particle velocity before the collision and bb `"v"_S` is the constant velocity of the sphere. Equation 14 is just an expression that the after velocity of the particle has the same angle of reflection with respect to the normal as its incident angle but the before velocity is modified by the sphere velocity. Of course, since the sphere has constant velocity (as if it is externally driven), the energy of the particle will not be conserved. Instead, the particles colliding with the front of the moving sphere will gain energy and those colliding with its back it will lose energy. Since the sphere is moving in the vertical (z) direction, the drag force on the sphere is the sum of the changes of the z component of the momentum of the particles hitting it during a single iteration time frame:

`F_d= sum_(i=1)^n -2 m_p [bb hat "n"_(ri) cdot [bb "v"_(bi)-bb "v"_S)] (bb hat "n"_(ri) cdot bb hat "z") \ \ \ (15)`

where `m_p` is the mass of a particle and the sum is over the total number of collisions, n, during a single time frame.

Description of the Graphics and the Animation

The corrugated sphere is treated mathematically as three dimensional (3D). But graphically it is not attempted to show its 3D nature. The column containg the sphere and the gas particles is also 3D mathematically with a width, W, depth, D, and height, H but I make no attempt to show it as 3D. Parallel to the gas column, I provide histograms of the particle density and particle energy Vs height, z. You should note that the density and energy tend to increase ahead of the sphere and decrease after the sphere. This is analogous to the bow wave and the wake wave that a boat creates. In addition, I provide time plots of the force on the sphere and the coefficient of drag computed from equation 2a. Both of these plots are time-lagged and averaged by 20 frames in order to average out large fluctuations or "volatility" as it's called in the stock market. Also the same data are printed out on Canvas 2.

Critique of Results

It is obvious from these results that the digital gas model that I have used is not adequate to see the expected values for the drag coefficient of the sphere. This hard sphere scattering gas model is perfectly adequate to yield the Maxwell-Boltzmann energy distribution. So what is the missing feature of the model? I think that an important clue comes from the fact that the hard sphere model also does not support propagation of a sound wave. Even with a repulsive potential instead of the hard sphere collisions, a propagating sound wave is not realized.

In order to get sound wave propagation in a digital model one must have both repulsive and attractive potentials between the particles. With this model there will always be an equilibrium distance between particles. Think of how we get a sound wave on a guitar string. The sound wave relies on the equilibrium tension of the string. If the tension is zero, then we can't induce a sound wave. Providing both repulsive and attractive forces between particles is similar to applying tension on a guitar string. It provides a potential derived built-in "pressure" which the hard sphere scattering model cannot provide. When we use a transducer at one end of a tube to induce a wave, the built-in pressure is modified because the distance between particles is is slightly changed with respect to equilibrium distance. This is equivalent to changing the displacement of a stretched spring which will always result in sustained vibrations in the absence of damping. Damping does not conserve energy and results in reduced number of vibration cycles. But the repulsive-attractive potential does conserve energy and therefore the sound vibrations will propagate.

Canvases 1 and 2: Animation of Sphere and Particles and Hisograms of Density and Energy
Canvas 3: Plots of Results as Well as Printed Data
Canvas 4: Plots Vs Sphere Speed of Average Force and Average Coefficient of Drag