Here I will examine the dynamics of a gas in a toroid where the gas has learner selectable initial rotational speed and density modulation frequency along the median (equator) of the toroid.
A toroid is an axisymmetric object that is characterized by two
parameters:
1. A median circle of radius `R_s` (similar to the equator of a sphere).
2. A set of tubular circles of radius `r_s` centered at this median .
To locate any point on the tubular circle one also needs to
define the azimuth, `theta`, pointing to a point on the median circle
as well as the altitude `phi` of the point on the tube.
Thus, an equation in Cartesian coordinates that defines the point can be
written
`bb "r"_v=[r_x,r_y,r_z]=[ R_s cos theta + r_s cos theta cos phi,\ \ \ R_s sin theta + r_s sin theta cos phi,\ \ \ r_s sin phi] \ \ \ \ \ (0)`
For viewing on the canvas, the toroid can be rotated by multiplying each toroid vertex position vector by a 3`x`3 rotation matrix, M. Particles are added to the toroid in groups at intervals along the `theta` angles. These Particles can be given a combination of random velocities as well equal drift speeds. For this animation. the drift speeds will be along the median or `theta` direction (i.e. they will be locally parallel to the toroid median line). For display only, the particles also need to be rotated with M in order to remain inside the toroid.
The initial velocity vectors of the particles are defined by the rotation speed `v_theta` and the random speed `v_r`.
`bb "v"=v_theta hat bb theta+(v_(rx)hat bb "x"+v_(ry)hat bb "y"+v_(rz)hat bb "z")`
where`hat bb theta=-sin theta hat bb "x"+cos theta hat bb "y"`
and`sqrt(v_(rx)^2+v_(ry)^2+v_(rz)^2)=v_r`
and the relative values of the three components are chosen using the random number generator: To start, define:`bb "v"_t=[v_(xt),v_(yt),v_(zt)]=2[rnd()-1/2,rnd()-1/2,rnd()-1/2]`
where each `rnd()` yields a different random number between 0 and 1. Now divide `bb "v"_t` by its magnitude to make it a unit vector:`hat bb "v"_t=bb "v"_t/|bb "v"_t|`
finally we have the vector version of `v_r`:`bb "v"_r=v_r hat bb "v"_t`
The Moving Gas inside the toroid can be thought of as a wind tunnel similar to what's used to test the lift and drag of aircraft bodies. Most wind tunnels have capabilities for very fast air speeds, `v_t`. The ratio of `v_t` to the speed of sound in still (non-moving) air is called the Mach number, M. The equation for the speed of sound in an ideal gas is usually given as:
`v_(sound)=sqrt(gamma (k T_K)/m)`
Here `gamma=C_P/C_V` is the ratio of the heat capacity at constant pressure to that at constant volume `k` is Boltzmann's constant, `T_K` is the Kelvin (absolute) temperature, and m is the mass of the particle. We have seen for all thermal problems that `k T_K` is exactly the same as the average random kinetic energy `k T_K=E_(avera"ge")`. So the equation becomes:`v_(sound)=sqrt(gamma (E_(avera"ge"))/m)`
where m is the particle mass, taken to be 1 in this program. Now the `random` kinetic energy is the same as the thermal kinetic energy in still air i.e. it does not include the rotational energy. The running rotational energy is a bit smaller than the starting rotational energy because the average particle distance from center increases due to the need to make circles in the toroid. Our toroid Mach number, M, can be written as`M=v_theta/sqrt(gamma (v_r^2)/2)`
For our simple spherical particles `gamma=5/3`.
Hereafter, a bold symbol denotes a vector which has [x,y,z] components. The particle's position can be characterized by two vectors: Its azimuth `bb "R"_theta` and its altitude `bb "r"_phi`. The actual position vector will be called `bb "r"_v` and it will be the sum of the two vectors:
`bb "r"_v=bb "R"_theta+bb "r"_phi \ \ \ \ \ (1)`
To determine whether the particle is close enough to the inner wall to be reflected, we need the magnitude of `bb "r"_phi`. We therefore solve equation 1 for `bb "r"_phi` and obtain:`bb "r"_phi=bb "r"_v-bb "R"_theta \ \ \ \ (2)`
Now we need the value of `bb "R"_theta`. We know that this vector neets the toroid median line at some angle `theta`. The vector `bb "r"` is defined by both a `theta` and a `phi` as can be seen by examining equation 0. Note that .`sin phi` =`r_z/|(bb "r")| \ \ \ \ (3)`
From equation 3 we can easily obtain `cos phi=sqrt(1-sin^2 phi)`. Now we can use equation 0 and `cos phi` to obtain `bb "R"_theta`.`bb "R"_theta=R/|bb "r"|[r_x/(cos phi),r_y/cos phi,0] \ \ \ \ \ (4)`
From equations 2 and 4 one obtains `bb "r"_phi` and this vector is along the unit normal `hat bb "n"` to the inside of the toroid.`hat bb "n"=bb "r"_phi/(|bb "r"_phi|) \ \ \ \ \ (5)`
When the particle collides with the inside, the velocity increment:`delta bb "v"=2 (hat bb "n" cdot bb "v") hat bb "n" \ \ \ \ \ (5)`
is added to the incident velocity in order to conserve angular momentum and energy.To get a better view of the toroid's angular rotation angle sliders are provided. The nominal axes are the x axis perpendicular to the screen, the z axis vertical and the y axis horizontal. These are the axis directions generally used in physics texts. When all the slider settings are zero, the toroid axis is initially along screen's vertical direction. The conserved quantities are the angular momentum and usually the kinetic energy. These are listed on Canvas 2. However, when a radial centripetal force is applied, kinetic energy is NOT conserved since radial potential energy is also changing then. There is a slight variation in the total angular momentum but that is just due to the limited numerical effort to integrate the particle velocities and positions. These can be easily improved by using a 4th order Runge Kutta integration method.
The particles' average distance from center oscillates after startup. The variation of this average distance as well as the average rotational speed are plotted on the bottom part of Canvas 2. Note that the oscillation frequency is a function of the rotational speed of the particles. This is reasonable because the centrifugal force needed to keep the particles inside the toroid is also a function of this speed. So what happens is the gas density increases as the distance from the tube center. Gas density increase causes the pressure to increase likewise. So the increased pressure acts like a spring and the oscillation frequency is a natural function of this spring. When we add more Radial Force Factor, these oscillation amplitudes are decreased.
The main feature of this animation is that the initial density of the gas is modulated versus angle `theta`. What is interesting is how quickly this modulation "washes out". The main reasons for reduction of the modulation are:
1. The adjustable thermal speed of the particles.
2. The collisions of the particles with mainly the outer inside boundaries of the toroid.
These collisions can be greatly reduced by applying a negative force (centripetal acceleration) along the radial vector of the particle. The magnitude of this force can be adjusted by changing the Radial Force Factor slider. With a value of 1, this parameter is exactly the force needed to keep the particles going in a circle at their initial rotational speed.The number of modulation cycles `m` Vs angle 'theta' can be adjusted using the Modulation Cycles slider. The depth of density modulation is not adjustable and is 2/3, meaning that the density varies from 1 to 1/3 of its peak value. On each iteration of the particle movement the number of particles Vs angle `theta` are placed in bins. The number of angle bins is 64 which is `2^6`. Some number `2^n` where n is an integer is the required number for the Fast Fourier Transform. (FFT) that is used here. The exponent n=6 is a convenient exponent for our bins here. The FFT works with both real and imaginary parts. In this case the real part input is the 64 bins particle numbers and the imaginary part input is an array of 64 zeros. The output of the FFT is an array of the frequency and phase components of the angle bin distribution. The largest element are the zeroth and first harmonics `A[f_0]` and `A[f_1]` where `f_0=0` and `f_1` is the number of modulation cycles. The phase of the first harmonic is computed using the equation `phi_1=arctan((A[f_1]_(imag))/(A[f_1]_(real)))` Then a smooth sinusoidal curve using the magnitude ratio ` |A[f_1]|/|A[f_0]|` and phase is plotted.
`f(theta)=a_0 |A[f_1]|/|A[f_0]|cos m (theta+phi_1)`
where a_0 is the initial modulation amplitude. The amplitude of the wave is proportional to the ratio of the magnitudes of the FFT first harmonic to that of the zeroth harmonic. When the animation is running, note that the curve amplitude and position changes to agree with the bin plot data.