One of the frequent topics of discussion on the larryniven-l science fiction email list is Larry Niven's novel, Ringworld, about a huge artificial habitat in the form of a ring about two astronomical units in diameter. It is heated by a central star, and spins for artificial gravity, with "spill mountains" that hold the air in. On this ring, made of a fantastic material called "scrith" (carrying astronomical stresses due to the artificial gravity), is a puncture called "the Fist of God" that was made by a meteor. Fortunately, the scrith is ductile enough for the Fist of God to have high rims to keep the air leakage rate down. This impact has led to discussions of how the Ringworld would respond dynamically to such an event. In the words of the old man in Harry Chapin's song, 30,000 lb. of Bananas, "Boy, it sure must've been somethin'." Obviously there would be tremendous shock waves, traveling at very high speed through the fantastically strong (and also fantastically stiff) scrith. But I don't know much about shock waves, and don't have material properties for scrith, and so I was wondering about the longer term response, the low frequency vibrations that would be excited by such an impact, and felt throughout the whole ring for a long time (and which I will argue are less sensitive to material properties).
To my structural engineer's eye, the Ringworld looks kind of like a violin string, except it is curved and has no ends. Shock waves (and other high speed waves) would involve axial stretching or shrinking of segments along the ring, like sound moving along a string between two Dixie cups, but the low speed waves would be lateral waves, like a violin string's response to being plucked. I want to know the shapes of these low frequency waves and their associated frequencies. This is called a modal analysis. I am interested in "natural" frequencies, which is how it vibrates after you let go of it, as opposed to "forced" vibrations, which are driven by some continuing external energy source.
Most of the vibrating systems I studied in school involved "simple harmonic motion," like a violin string or a pendulum on the equator. This is motion that alternates between zero velocity at maximum displacement and zero displacement at maximum velocity. The motion usually takes the form of a sine wave, varying in time in proportion to the sine of (the frequency times time plus an offset). The offset, or phase angle, lets you match boundary conditions on both initial position and initial velocity. If there is damping, an energy dissipation mechanism like the shock absorbers in a car, the motion will decay exponentially as well as sinusoidally, and if there enough damping ("critical" damping), it will cease to be sinusoidal at all.
All of the systems I'm going to be describing are linear in the sense that the forces acting on each mass are linearly proportional to the various displacements, velocities, or accelerations. If there is a non-linear effect, such as the fact that the lateral motion of a violin string depends on how tightly it is stretched axially, I will linearize the system by assuming that there are only small motions about an equilibrium state (the violin string is already stretched before I start calculating a small lateral motion). Otherwise, the problem is too hard for an amateur like me to deal with. Also, because the Ringworld is rotating, in order to keep the problem simple enough for me to deal with, I need to pull some math tricks that I can't pull if I consider damping. So henceforth I will ignore damping (which for aerospace structures engineers is usually an afterthought anyway).
In order to describe the vibrations in terms of small motions about an equilibrium, I pretty much have to use a coordinate system that rotates with the Ringworld. This forces me to consider "ficticious" forces such as Coriolis force and "centrifugal" force; in reality, objects in a centrifuge are accelerating inward, but in a rotating coordinate system, they may appear stationary. Generally speaking, the equation that relates accelerations in two different coordinate systems has five terms. If I have a point, P, an "accelerated" reference frame, B, and an "unaccelerated" reference frame, A (the kind I want), the acceleration in frame A includes:
- The acceleration of point P in frame B.
- The acceleration of the origin of frame B in frame A (zero in this case).
- The effects of changes in the rotation rate of frame B (zero in this case).
- Centrifugal force: omega cross (omega cross r) for those of you who are into cross products. "omega" is the rotation rate of frame B, and "r" is the position of point P in frame B.
- Coriolis force: 2 * omega cross velocity (in frame B). This is the least obvious of the acceleration terms, and the easiest one for engineers to overlook, which is why M. Coriolis is so honored. In terms of another Niven novel, The Integral Trees, "Forward takes you out, out takes you back, back takes you in, and in takes you forward." Another way of looking at this is that if you are moving radially outward in a rotating centrifuge, your true velocity is increasing (in the spinward direction), and this requires a force, even if your velocity appears constant in the rotating reference frame. If you are moving in the tangential direction (spinward), your velocity also contributes to the centrifugal effect (which is why the cross product formula for Coriolis force still works).
So we're in a rotating coordinate system. Forget the violin string analogy for a moment and imagine that you're riding a merry-go-round while trying to play tetherball. Let's also imagine that the merry-go-round has curtains around the perimeter so stationary objects outside it don't give us another reference frame.
If you pull the tetherball outward from its equilibrium position (away from the merry-go-round's center) and push it forward (in the direction the merry-go-round is moving), the tetherball will appear to move in a circle. This is shown in the figure,
The large dotted circle shows the tetherball's equilibrium point moving relative to the merry-go-round's pivot. The long red line connects the pivot to the equilibrium point. The short red line connects the equilibrium point to the ball's current position, the arrow shows the ball's current velocity, and the small circle shows the ball's path around its equilibrium point. As is true of all of these figures, the displacements have been scaled to make them look good, so don't try to read anything into the scale; the equations have been linearized, so the scale is fundamentally arbitrary.
As the ball swings, its momentum (the centrifugal force of the ball moving relative to its mount) holds it away from its equilibrium point as the restoring force in the tether pulls it inward. Note that this motion is not harmonic; the ball does not pass through its equilibrium position, but circles around it. If we try to force it to pass through its equilibrium position while in motion, something more complicated will happen which I will get to shortly. If I want to write an equation for the position of the ball at any time in its cycle, I need more than one displacement vector and one sine wave. I will want another displacement vector from when the ball is a quarter of the way around the circle, and another sine wave that is 90 degrees out of phase with the first. Then I will be able to express the ball's position as the sum of two fixed vectors and two sine waves. In some ways it is easier to think of these two vectors as different phases of the same mode of vibration, but in other ways it is easier to think of them as a pair of different modes that happen to share the same frequency. I show this in my next figure,
We could just as easily push the ball backwards, opposite to the merry-go-round's rotation. But there is a difference. In the forward case, the ball's rotation and the merry-go-round's rotation add. In the backward case, they subtract. So the periods of the motions in the two directions will be different. In the backwards case, the ball needs more speed in order to fight the tether's restoring forces, and so it will oscillate with a higher frequency (the motion will have a lower period). For the simple example problem I solved which trying to debug my Ringworld analysis program, the merry-go-round spun at 2 radians per second, the ball weighed 1 kilogram, and the springs that held it in place (in two perpendicular directions) had stiffnesses of 104 Newtons/meter. The frequency of motion in the posigrade direction already shown was 8.198 radians/second and in the retrograde direction was 12.198 radians/second. The second pair of vectors is shown in the next figures,
In real life, we usually don't see pure modes like these circles, but combinations of several modes at once. You'll notice that to get these pure modes, I had to displace the tetherball and then push it forward or back. If I had just pushed it without first displacing it, I would have excited two modes at once, and would have gotten a more complicated motion that did swing back through the equilibrium point, but didn't quite repeat itself on the next swing. The two modes vibrate at different frequencies, and so the tendencies to swing clockwise and counterclockwise don't quite cancel. We then see, in an exaggerated form, the familiar motion of a Foucault pendulum mounted well away from the Earth's equator:
Now let's talk about the Ringworld. I will assume that the radius is exactly 1 astronomical unit (AU), and the central star's mass is exactly that of our sun. I will also assume that the artificial gravity level is exactly 1 standard gravity. I will assume the total mass of the Ringworld is the same as that of the planet Jupiter, and that the long-period vibrations that I care about involve motion that is entirely within the plane of the Ring (It is over a million miles wide, but only a few meters thick, so I think this is a safe assumption). Some more assumptions may sound a bit strange: because the Ring is so thin relative to its diameter, I will take the violin string analogy very seriously, and assume that the Ring's bending stiffness within its own plane is negligible. That is to say, the Ring's resistance to lateral displacement is due entirely to the stiffening effect of the astronomical axial load it is carrying, and the stiffness of scrith (Young's modulus) is unimportant. But there are also higher-frequency vibration modes involving stretching and shrinking of segments of the Ring, so I need to assume something about the axial stiffness. Specifically, despite having just said that the in-plane bending stiffness is so low that it is unimportant, I will now assume that the axial stiffness is so high that it is unimportant. That is, the "high-frequency" modes are so high that they don't relate in any important way to the low-frequency modes that I'm trying to calculate. In one direction, I set the stiffness to zero (apart from the stiffening effect of the high tension load), and in the other direction I set the flexibility to near zero (stiffness to near infinity). You can think of the Ringworld as a musical instrument: I'm trying to predict what the bass will sound like, and I don't care if the treble is wrong. By the way, the stiffening effect of preload is proportional to the Ringworld's mass, and the frequencies are functions of the ratio between the two, so the actual mass I assume isn't important, either. To do a modal analysis using a "finite element model," I have to break the ring into discrete pieces, the nodes at which I measure movement and the structural elements that connect them. As of this writing, I haven't played with models that had more than 24 nodes, but for the lowest vibration modes, which I am interested in, this seems to be enough. You can think of my model of the Ringworld as a well-oiled bicycle chain with 24 links in it.
A warning about the Ringworld: it is gravitationally unstable, and requires a control system that had to be repaired in a sequel called The Ringworld Engineers. (Hey, Ted or Andy, isn't there a page on Known Space that explains this?) Even if it were neutrally stable, as a free-floating object, it would have "rigid-body modes" that have frequencies of zero and involve no additional stresses beyond equilibrium stresses. In three dimensions, there are six rigid-body modes: three uniform motions (translations) and three rotations. In two dimensions, there are two translations and one rotation. For my purposes, these are merely a nuisance.
Besides the gravitational instability, have I introduced an additional instability by zeroing out the bending stiffness? In other words, if I spin up a bicycle chain in free fall, will it form a circle, or will it stretch into a doubled-up straight line (I'm too lazy to consider any intermediate possibilities)? Mother Nature wants to minimize energy (0.5 * spin rate * moment of inertia squared) while conserving angular momentum (spin rate * moment of inertia). A quick check of the moment of inertia of a circle (mass * radius squared) and a line (mass * 1/12 * length squared), where length = 0.5 * circumference, shows that the circle is more stable, so my assumptions have passed at least one reality check.
Now let's think about that centrifugal force (omega cross (omega cross r)) and the astronomical forces in the scrith resisting it. If the mass per unit arc length is rho, the equilibrium hoop load will be rho * (r * omega) squared (the corresponding "preload" parameter in my finite element is a hair different because my model is discrete rather than continuous). Scrith has to have a fantastically high strength to weight ratio, but are there also limits on its possible stiffness? Let's assume the Ring has a given angular momentum, and the control system does not add angular momentum to try to maintain a fixed rotation rate. In the limiting case of very low stiffness (but no limit on how far it can stretch), the Ring becomes a large bunch of nearly independent pieces flying off on trajectories that are nearly hyperbolic orbits. (It has to be rotating much faster than escape velocity in order to have 1 G of artificial gravity, 1211 km/sec as opposed to 29.8 km/sec. Actually, the official number from the novel is 770 miles/second, or 1239 km/sec. So shoot me.) But this kinetic energy is finite, and if there is no limit on how much strain energy the scrith can store, there is no problem. Eventually the restoring forces in the scrith will overcome the Ring's outward momentum. (If the radius is varied as mass and angular momentum are held constant, the hoop load is inversely proportional to the radius cubed.)
On the other hand, if the control system maintains a constant angular rate as the scrith stretches, we can compare the equation for the hoop load as a function of radius, and the equation for strain (fractional change in circumference) of the scrith as a function of hoop load and stiffness. If the radius of the Ring before it was spun up for artificial gravity was R0, the cross section area was A, and Young's modulus (stiffness) for scrith is E, then the equilibrium radius of the Ring becomes
r = 2*Pi*E*A*R0 / (2*Pi*E*A - mass * omega squared * R0) .
This blows up as E*A drops toward mass * omega squared * R0 / (2 * Pi).
We could also look at what happens if the control system holds the artificial gravity level constant rather than the angular rate. I get
r = R0 * (1 + mass * perceived gravity / (2 * Pi * E * A) ),
so this would not bound the stiffness.
Since I am trying to linearize the Ringworld's motion in a reference frame that rotates at a constant rate, I am more interested in the calculation for constant rotation rate. I have already argued that physically, there is no problem with the scrith being flexible. Even if the control system tried to maintain a particular angular rate, it seems unlikely that it could respond fast enough (or be stupid enough) to exacerbate a vibration mode involving the overall stretching of the Ring. Do I care about this supposed lower limit on scrith stiffness? Actually, I do. I am incorporating tidal effects (the reduction in gravity as a piece of the Ring deflects outward) in my system of equations for the Ringworld's motion, and also the increase in centrifugal force with increasing distance from the center. I already have three rigid body modes (two of which are downright unstable), and now the apparent tidal forces (including centrifugal force variations) give the appearance of being destabilizing in yet another way. I think I need my axial stiffness to be high in order to avoid numerical problems--but if it's too high relative to the lateral stiffness (due to preload), I will have a different set of numerical problems.
I tried to estimate the first (lowest) natural frequency by hand, ignoring Coriolis forces, and assuming that that the first mode shape was an ellipse, and that the behavior of a small area near the major axis of the ellipse was typical of the entire ellipse. I got an answer of 5.2 days for the period. Actually, I hadn't realized that the motion wouldn't be harmonic, and so my guess was that this mode would be an ellipse with the major and minor axes always pointing in the same two directions (say, 0 and 90 degrees), with the deflections passing through zero as the two axes reversed. I also expected another mode at the same frequency with the axes offset by 45 degrees (linear combinations of these modes allowing any desired orientation).
I also tried to do a simulation using a fourth-order Runge-Kutta differential equations solver, but I concluded that my equations were "stiff" in the numerical sense of being badly behaved, too sensitive to error to give me useful results. This is what I meant about the axial stiffness being too high relative to the lateral stiffness. Rather than trying to use a better diff. eq. solver, or a better way of defining the variables in the problem, I gave up and decided to do a proper modal analysis instead. (I also thought later about using a Rayleigh-Ritz approach (Fourier series), but by then I'd invested a lot of time in the finite element program, and wasn't sure anything else would work any better.)
It soon dawned on my that I had studied single-degree of freedom mechanical vibration problems in detail in school, and had studied matrix methods of static structural analysis in detail, but that I really didn't understand how to work dynamics problems that were this complicated. So I ambushed a structural dynamicist (who shall remain anonymous) at my favorite coffee bar and borrowed a dynamics textbook (by Leonard Meirovich).
The equations of motion in "problem space" are
M * x-dot-dot + G * x-dot + K * x = 0, where M = mass, G = Coriolis, and K = stiffness matrices, and x is the displacement vector.
Meirovich transforms these from a second-order system of differential equations (problem space) to a first order system ("state space") by defining a new X vector twice as big, with both displacement and velocity terms, and combining K with the other matrices to get a symmetric matrix, M*, and a skew-symmetric matric, G*. Then comes a bunch of voodoo involving a Cholesky decomposition (expressing the M* matrix as the product of a triangular matrix and its transpose), some linear algebra with complex numbers, and a change of variables. Eventually, this turns into an eigenvalue and eigenvector problem for a matrix, A, that is real, symmetric, and positive-definite, which was fine with me because I had a subroutine that I copied out of a numerical methods book years ago that could handle this class of problems (Jacobi's method). Meirovich goes to a lot of trouble to be able to work with real, symmetric matrices (so the eigenvectors are real), but the voodoo requires G to be skew-symmetric. A symmetric damping matrix, C, would go here if this were a problem in which damping is important, but that requires an incompatible brand of voodoo, so we won't consider damping. With no damping, the solutions come out like the tetherball problem: sane people talk about vectors of real numbers times sines and cosines of frequencies times time, and mathematicians talk about pairs of complex conjugates times e to imaginary powers. But either way you look at them, the solutions have to come in pairs that have the same frequency, and all of these vectors should be orthogonal (given appropriate transformations).
The Meirovich voodoo requires M and K to be positive-definite. I tried several shortcuts to getting K (the stiffness matrix) to be positive-definite, adding stiffness terms to represent various kinds of springs and mechanisms. Still I was getting erratic results: frequencies that didn't come in pairs (a typical sign of an ill-conditioned matrix), and mode shapes that looked bizarre. I tried using a gradient routine to improve on the accuracy of the Jacobi routine, but it didn't help. Next I used a formal approach to eliminating the rigid-body modes: writing constraint equations (various linear combinations of motion were set to zero), dividing my problem space variables into "independent" and "dependent" degrees of freedom, and eliminating the dependent ones (transforming my M, G, and K matrices into the constrained space). (This was easier than it sounds because I had an old homework-related program that did exactly what I wanted.) That helped somewhat. Next I tried to get rid of the excessively stiff axial motions by adding more constraint equations, representing rigid links in my "bicycle chain" rather than very stiff ones (and thus eliminating these now-dependent degrees of freedom from the problem space). That helped a lot, but still the results were somewhat erratic, especially when I started playing with the number of nodes in the model. Odd numbers of nodes gave especially bad results. I had another old homework-related routine that verified that my A matrix was positive-definite (indicating all positive eigenvalues), but for a 21-node model the Jacobi routine would return a negative first eigenvalue (which sucks, because the frequencies of vibration are the square roots of the eigenvalues).
My next move, an act of desperation, was to switch from treating mass as being "lumped" at the nodes to being distributed along the elements ("coupled mass"). This actually helped noticably. One wrinkle was that I forgot to apply the "coupled" scheme to the Coriolis matrix, then discovered and fixed the error, but couldn't find any difference in the output. It turned out that the effect that this coupling had on the Coriolis matrix was removed as a side effect of applying the rigid bar constraints.
Finally I decided to see if it made any difference whether I used Jacobi's method or the inverse power method, and the first frequencies indeed were much better behaved. I still use Jacobi's method, but now I invert my "A" matrix first (and then invert the eigenvalues). This way, whatever problems there are with ill-conditioning manifest themselves in the high-frequency modes that I don't care about.