Equilibrium figures of spinning bodies with self-gravity

Equilibrium figures of spinning bodies with self-gravity

Icarus 172 (2004) 272–303 www.elsevier.com/locate/icarus Equilibrium figures of spinning bodies with self-gravity Keith A. Holsapple ∗ Department of ...

1MB Sizes 1 Downloads 73 Views

Icarus 172 (2004) 272–303 www.elsevier.com/locate/icarus

Equilibrium figures of spinning bodies with self-gravity Keith A. Holsapple ∗ Department of Aeronautics and Astronautics, University of Washington 352400, Seattle, WA 98195, USA Received 13 February 2004; revised 28 May 2004 Available online 3 September 2004

Abstract The study of the equilibrium and stability of spinning ellipsoidal fluid bodies with gravity began with Newton in 1687, and continues to the present day. However, no smaller bodies of the Solar System are fluid. Here I model those bodies as elastic–plastic solids using a cohesionless Mohr–Coulomb yield envelope characterized by an angle of friction. This study began in Holsapple 2001. Here new closed-form algebraic formulas for the spin limits of ellipsoidal shapes are derived using an energy method. The fluid results of Maclaurin and Jacobi are again recovered as special cases. I then consider the stability of those equilibrium states. For elastic–plastic solids the common methods cannot be used, because the constitutive equations lack sufficient smoothness at the limiting plastic states. Therefore, I propose and study a new measure of the stability of dynamic processes in general bodies. An energy-based approach is introduced which is shown to include stability approaches used in the statics of nonlinear elastic and elastic–plastic bodies, spectral definitions and the Liapunov methods used for finitedimensional dynamical systems. The method is applied to spinning, solid, strained bodies. In contrast to the special fluid case, it is found that the strain energy term of solid materials generally induces stability of all equilibrium shapes, except for two possible cases. First, strain softening in the elastic–plastic law can result in instability at the plastic limit spin. Second, a loss of shear stiffness can give unstable states at specific spins less than the limit equilibrium spins. In the latter case, a solid spinning ellipsoidal body without elastic shear stiffness can spin no faster than with a period of about 3.7 hr, else it will fail by shearing deformations. That is distinctly slower than the oft-quoted limit of 2.1 hr at which material would be flung off the equator by tensile forces. However, the final conclusion is that neither cohesion nor tensile strength is required for the shapes and spins of almost all of the larger observed asteroids: we cannot rule out rubble-pile structures.  2004 Elsevier Inc. All rights reserved. Keywords: Asteroids; Asteroids, rotation; Comets; Equilibrium shapes; Equilibrium figures; Stability

1. Introduction Possible shapes and spins of the bodies in the universe are restricted by their composition, and by their gravitational forces. There is a classical literature on equilibrium shapes for inviscid incompressible ellipsoidal fluid bodies, with the primary application being to stars. The extensive monograph “Ellipsoidal Figures of Equilibrium” by Chandrasekhar (1969), presents the early contributions of Newton, Maclaurin, Jacobi, Meyer, and Louiville and others mentioned below; the more recent review Lebovitz (1998) gives some recent developments. The interested reader can

refer to those references for background.1 Newton considered only slightly oblate fluid bodies. Maclaurin discovered combinations of spin and arbitrary oblateness that were in equilibrium; those are now called “Maclaurin spheroids.” Jacobi discovered certain equilibrium ellipsoids with axes of unequal length; those are now called the “Jacobi ellipsoids.” Each Maclaurin or Jacobi body has a single possible equilibrium spin. Then internal motions in fluid bodies relative to a rigid rotating frame were introduced into the theory. Dirichlet analyzed ellipsoidal bodies with motions that were linear functions of position in a rotating Cartesian coordinate system. After the death of Dirichlet, Dedekind published the results

* Fax: 206-543-0217.

E-mail address: [email protected]. 0019-1035/$ – see front matter  2004 Elsevier Inc. All rights reserved. doi:10.1016/j.icarus.2004.05.023

1 The specific references to the older literature are not reproduced here.

Equilibrium figures of spinning solid bodies

of Dirichlet with additions. Riemann showed that such allowable motions include rigid spins and uniform internal vorticity. That family of shapes and motions are now called “Riemann ellipsoids.” Thompson (Lord Kelvin) and Tait suggested that these fluid bodies might evolve into a nonellipsoidal shape and then into a planet–satellite system or into a binary star. Poincaré discovered a sequence of “pear-shaped” equilibrium configurations, and hypothesized stable evolution into planet–satellite systems. Darwin suggested evolution into binary star systems. This “fission theory” was studied by Darwin, Liapunov and Jeans for a number of years, but came to a halt in 1924 when Cartan showed that the pear shapes were unstable. Roche considered the stability of a small satellite rotating around a rigid spherical planet, leading to the famous “Roche limit.” Studies about fluid bodies continue to the present. Lebovitz (1998) reviews the recent work. Fasso and Lewis (2001) present a detailed numerical study of the Riemann ellipsoids, and find that there are complicated regions in shape space where the Riemann ellipsoids are spectrally stable, and further, that all satisfy a more general stability criteria over long time scales. The main interest here is for Solar System bodies. With the exception of the four gas giants, the bodies in the Solar System are solid: they are composed of rocks, dirt, ices, or metals. The fluid studies can give only limited insight into the equilibrium shapes and stability of asteroids and comets, since the assumption of fluidity is entirely inappropriate. Studies using material models that are appropriate for rocks, soils and other solid materials are required. The analysis of solid behavior introduces significant complexities into the problem. A complete stress tensor is needed. The analysis of deformation needs both spatial and reference configurations. Stability considerations involve second-order terms in constitutive equations, so large strain and deformation constructs are needed.2 Failure criteria for general three-dimensional stress states are required for these complex materials. Fortunately, one can at least dispel with internal motions, and consider only limit states with rigid spins.3 I presented a study for solid bodies in (Holsapple, 2001), where I derived limits on equilibrium configurations for a cohesionless4 solid ellipsoidal body. The analysis determined the limiting internal stress distributions due to the forces of 2 Since theories of general continuum mechanics, plasticity and soil mechanics may not be well known for some readers of this paper, some background will be presented here. 3 (Lebovitz, 1998) explains that, in the 1880’s and for over 40 years, the solutions with internal motions were largely ignored, probably because of the prevalent belief that stars had sufficient viscosity to enforce rigid motion. That belief was later dispelled, so the internal motion cases again became relevant. Here, coming full circle, the stiffness of a solid body can be invoked to make those cases irrelevant. 4 One with zero tensile strength. Note however, that does not mean strengthless, see Section 4.

273

self-gravity and spin. The bodies were modeled as elastic– plastic materials with a Mohr–Coulomb yield model characterized by its angle of friction. That model is appropriate for “rubble-pile” structures that consist of a continuum assemblage of particles held together at moderate densities by mutual gravitational forces. It is now suspected that such a characterization includes most asteroids and perhaps comets. A major observation supporting that rubble pile interpretation is the fact that all larger asteroids have a spin within certain theoretical speed limits. A “rotation rate barrier” for rubble pile asteroids is often stated as about 2.1 hr/rev (Harris, 1996); that estimate is the result of a simple comparison of the gravity forces and the centripetal forces at the equator of a spherical body. The (Holsapple, 2001) results show that those “speed limits” for a spherical body are too high, and that the speed limits are constrained by the shear stresses, not the tension at the equator. Failure by shearing happens before material can be flung off of the equator. Further, the theory shows how the spin limits decrease in a definite way for elongated bodies. Here a new approach and substantial additions to that previous theory of (Holsapple, 2001) are presented. I use an energy formulation for the study of equilibrium for general material behavior. Using that formulation, new closed-form algebraic expressions are derived for the limit equilibrium states of cohesionless bodies, with explicit algebraic dependences on the ellipsoidal shape and the friction angle. The constitutive model for stresses within yield is basically irrelevant, so can include any nonlinearity. Since the limiting case of zero friction angle includes the fluid model, the classical fluid equilibrium states and formulas for the Maclaurin spheroids and the Jacobi ellipsoids are recovered as special cases in a new way. Thus, the results here are a major generalization of the fluid equilibrium results to include much more general materials including both fluids and solids. For fluids, it is known that many of the equilibrium states are not also stable, so are not possible physical states. Bifurcations in shape are also possible at certain spin rates, the primary one being that from the Maclaurin to the Jacobi bodies. Let α = c/a denote the ratio of the ellipsoidal shortest semi-axis c to the longest a. Chandrasekhar (1969) re-derives the result that the Maclaurin Spheroid states with α less than 0.3033 are not stable. Further, he finds that while cases with α between 0.3033 and 0.5827 are stable, that stability is lost if there is any amount of viscosity. Thus, many of the equilibrium states for a Maclaurin Spheroid are not stable. At the state with α = 0.5827, the bifurcation into a stable Jacobi ellipsoid is possible. For the Jacobi ellipsoids, (Chandrasekhar, 1969) finds that the equilibrium states are not stable if α is less than 0.3451. Since the solid equilibrium states include as limits the fluid states, it is important to consider the stability and bifurcations states of more general materials. When are they stable, and when are they not? What are the defining criteria for the loss of stability? For bifurcation? The fluid models used have zero friction angle, are incompressible, and

274

K.A. Holsapple / Icarus 172 (2004) 272–303

have zero shear rigidity. Which of those characteristics affects stability or allows bifurcation? What other aspects of the material model may influence stability limits? Are all solids stable in all limit states or not? Are bifurcations possible? Because of the differences between the allowable motions of a general elastic–plastic solid compared to a fluid, the classical, small oscillation perturbation methods to test stability in fluid analyses cannot be used. General stability criteria for a general material at some defined “equilibrium state” should investigate the boundedness of the motion resulting from all possible perturbations, but solutions are not usually available for those general cases. Therefore I propose a new definition of stability based on energy considerations. The definition includes known requirements for static stability for general solid materials, is sufficient to recover the bifurcation states and some stability limits for the fluid bodies, and shows what material behavior is required to give stability and bifurcation possibilities for solid bodies. It is found that the admissible range of spin rates may be substantially restricted by the stability considerations, but that depends significantly on the material properties; which are highly uncertain. For those stability limits, the associated spin rates are distinctly less, and the periods greater, that those possible considering equilibrium only. The “rotation rate barrier” for rubble pile asteroids, often stated as about 2.1 hr/rev (Harris, 1996) can decrease to a spin with about 3.7 hr/rev for spherical bodies, and the limit has an even greater period (slower) for elongated bodies. However, the new stability limits are still above the majority of data points for the asteroids. This new more general analysis does not rule out most asteroids as candidates for the rubble pile structure. An outline of the paper is as follows. Section 2 presents background material from continuum mechanics that may not be entirely familiar to some readers. It gives the notation, the equilibrium equations, the virtual work theorem, and a theorem for average stresses; and introduces the special cases of body forces from potential functions, of elastic materials, and of spinning bodies, in order to relate to more common but special forms. Section 3 develops a new and general approach and definition of the stability of equilibrium states for general bodies, using energy forms. That is necessary since the more usual approaches are insufficient for the present application. It then relates this new approach to more traditional but special common forms. Section 4 develops the theory of plastic bodies used for soil and rock mechanics, especially the Mohr–Coulomb approach. It presents some data for dry soils in order to illustrate the general features and typical numerical values for the behavior of dry granular materials. Section 5 then focuses on spinning ellipsoidal bodies, giving the body forces from spin and selfgravity, the kinematics of spin, the average stresses which are independent of the constitutive model, and the special forms for the virtual work, and stability definition and the stress work term for a solid. Then the derivation of the max-

imal spin states as a function of the material properties of a Mohr–Coulomb material and the ellipsoidal shape in Section 6. The major result is Eq. (6.3), a simple but general algebraic form for the limiting spin states. Section 7 considers simple fluid bodies, in order to show how the classical fluid Maclaurin and Jacobi figures are included as special cases, and derived here in a new way. Section 8 then moves on to the more general solid bodies and presents families of curves for limiting spin as a function of the shape and angle of friction. The data for the asteroids are compared to these limiting curves. It presents discussion of the stability of equilibrium for those solids and identifies certain material behavior that can lead to loss of stability for spin states near the equilibrium limits. Finally, Section 9 gives some concluding comments. Appendix A develops various algebraic forms for the stability criteria that are used in the paper.

2. Equilibrium, virtual work, average stresses, and stability Some notation and results for a general continuum will be presented in this section. Special results for a spinning body are derived. 2.1. Preliminaries and notation The standard notation of continuum mechanics will be used (Truesdell and Noll, 1965, hereafter referred to as T&N). A motion is the time dependent mapping of material particles with reference coordinates X to spatial positions x: x = χ(X, t).

(2.1)

The mapping χ at any given time is called a configuration of the body. The displacement vector is u = x − X.

(2.2)

The deformation gradient tensor is   F = Grad χ(X, t) ,

(2.3)

where Grad refers to derivatives with respect to the reference coordinates X. The determinant of F is denoted as J . Integrals over the reference configuration use a volume differential dV , those over the spatial configuration use dv, they are related by dv = J dV . The velocity and acceleration fields are given as v=

dχ(X, t) , dt

(2.4)

d 2 χ(X, t) (2.5) . dt 2 The Cauchy stress tensor is denoted by T (also below by σ ), and the first Piola Kirchoff stress tensor by TR . They are related by T = J −1 TR FT . Since the theory does not include

a=

Equilibrium figures of spinning solid bodies

couple stresses, the Cauchy stress tensor is symmetric, but the Piola Kirchoff tensor is not. I adopt the general constitutive assumption that, for each particle X, the history of the deformation gradient over all past time τ  t determines the stress tensor T at time t, this can be written as5   T(X, t) = T F(X, τ ) , τ  t, (2.6) where T is a functional of the past deformation gradient history. For some materials the domain of deformation gradient histories may be restricted in some way,6 and then there can arise a part of the stress not determined by Eq. (2.6). For specific results below I assume a particular plasticity model. The body force per unit mass is denoted as b and the mass density in the spatial configuration as ρ. Then Cauchy’s first law of motion is given as

where dm = ρ dv is the mass differential, T : δe denotes Tr[Tδe], the trace of the tensor product Tδe, t = Tn is the surface traction vector on the outer surface with outward unit normal n, and δe is the symmetric part of grad[δu]; i.e., the small strain tensor of the virtual displacement field δu. It is also given as the symmetric part of δFF−1 . The symmetry of the Cauchy stress tensor has been used. The surface work term is zero if the surface is not loaded, which is the case of main interest here. Equation (2.9) is entirely general and must hold for any smooth field δu. If the displacement increments are those along an actual history of the body, then they can be related to the present velocity field as δu = vδt. In that case Eq. (2.9) gives a result for time derivatives along the motion:    b · v dm + t · v da − T : e˙ dv Mass

div(T) + ρb = ρa,

(2.7)

where div refers to derivatives with respect to the spatial coordinates x. Alternatively, it is expressed in the reference configuration by using the Piola–Kirchoff stress and the mass density ρ0 in the reference configuration as   Div TR + ρ0 b = ρ0 a, (2.8) where Div refers to derivatives in the reference configuration. I shall refer to fields that satisfy Eqs. (2.7) or (2.8) as being “in dynamic equilibrium.”

275

Surface



Vol

d a · v dm = KE, dt

= Mass

where KE is the kinetic energy of the body, a form referred to as the theorem of virtual velocities. The calculation of the virtual work can also be made using Eq. (2.8) with the Piola stress and integrating over the reference configuration, the result is   b · δu dm + tR · δu dA VW = Mass

Surface



Mass

Surface



=

a · δu dm,

T : δH dV



2.2. Virtual work and energy theorems There are many work and energy results looking similar but which are actually different, and for different applications. I will present some of them in some detail to clarify the context, to pay special attention to terms arising from spins, and to compare to the related theory of finite-dimensional systems where they are well known. Taking the dot product of Cauchy’s first law, Eq. (2.7), with a small “virtual displacement”7 δu(x), integrating over the body, and using the divergence theorem, gives the general theorem of virtual work:    VW = b · δu dm + t · δu da − T : δe dv Vol

(2.9)

Mass

5 Such a material is called a “simple material” and is general enough to include all common linear and nonlinear theories of elasticity, plasticity and viscosity. See T&N. 6 For example, incompressible materials only admit isochoric motions, and then there is an arbitrary pressure not determined by the constitutive equation. Instead there is a constraint on allowable motions. 7 It will be important to recognize results that hold for any displacement increments (virtual) and those that only hold for “actual” ones. Actual ones must satisfy all equations of the theory.

(2.10)



R

T

Vol

=

a · δu dm,

(2.11)

Mass

where δH = Grad(δu) is the material gradient of the virtual displacement field and tR is the traction with respect to the reference configuration. Equation (2.11) can also be obtained term for term from Eq. (2.9) by using the transformations between T and TR , δH and δe, dv and dV , and the area differentials da and dA. Alternatively, since any reference can be used, if the present state is taken as the reference then the instantaneous Piola stress is just the Cauchy stress, and Eq. (2.11) is just a different notation for Eq. (2.9). It is often the case that certain terms in the virtual work can be expressed as a first variation in some scalar global measure of the body. Such cases rely on additional and special assumptions, some of which are now considered as examples. 2.2.1. Conservative body forces Often the body force is the gradient of a potential function U (x):   b = − grad U (x) . (2.12) There are two different cases of interest. If the potential function is from an external field and is not affected by a change

276

K.A. Holsapple / Icarus 172 (2004) 272–303

in configuration of the body, then the potential energy of the body is defined as  U (x) dm. PE = (2.13) Mass

However, for the case of self-gravity, any change in the configuration will change the potential itself. In this case, the potential energy is defined (Chandrasekhar, 1969, Eq. (12)) as  1 PE = (2.14) U (x) dm. 2 Mass

In either case, corresponding to any increment δu in displacement, the body force term in the virtual work expression of Eq. (2.9) becomes the first variation in the potential energy:8  b · δuρ dv = −δPE (2.15a) Vol

and, for an actual motion, the body force term in Eq. (2.10) becomes  d b · vρ dv = − PE. (2.15b) dt Vol

2.2.2. Hyperelastic materials Another common case is for a hyperelastic material (T&N, Eq. (82.10)), then the Piola stress is given in terms of a strain energy function Σ(F) ∂Σ(F) (2.16) . ∂F Then the stress term in the virtual work Eq. (2.11) is   ∂Σ R T : δHT dV T : δH dV = ρ0 ∂F Vol  = δ Σ dm = δSE, (2.17a) TR = ρ0

where the (total) strain   energy of the body is defined as SE = Σ dm. The term Vol T : δe dv of Eq. (2.9) is also δSE. For an actual motion the stress work term in Eq. (2.10) is   d ˙ T dV . SE = T : e˙ dv = TR : H (2.17b) dt Vol

Vol

Generally, materials are not hyperelastic. However, in all cases, for any δu(x) there is always a well-defined increment of stress work defined as   δSW = T : δe dv = TR : δHT dV (2.18)

8 It is always assumed that the variation in density is so that the incre-

ment of mass is zero: δ(ρ dv) = δ(dm) = 0.

although this increment it is not generally a first variation of a strain energy functional of only the present configuration.9 Since the stress may always be considered the sum of a hyperelastic part Te plus a remainder dissipative stress Td , the stress work increment can be also decomposed as a sum of an increment of elastic strain energy plus an increment of stress work of the dissipative stresses defined as:   δSE = Te : δe dv, (2.19) δSW d = Td : δe dv so that   T : δe dv = TR : δHT dV = δSW = δSE + δSW d , (2.20) which can be used in Eqs. (2.9) or (2.11). Also, along any actual motion,10   ˙ T dv = d SW = d SE + d SW d , T : e˙ dv = TR : H dt dt dt (2.21) which can be used in Eq. (2.10) to give, whenever there is a potential for the body forces,  d d [KE + PE + SE] = t · v da − SW d . (2.22) dt dt 2.2.3. A kinetic energy term for freely spinning bodies Equations (2.9) and (2.11) concern only the current condition and arbitrary displacement increments for a body. A natural question is whether the terms involving the acceleration can be written as a variation in some energy such as the kinetic energy. However, the kinetic energy involves dynamics and the velocity, so an increment in the kinetic energy would generally be determined by some virtual increment in the velocity: an entirely independent consideration from the virtual displacements. However, if increments in velocity are related to increments in displacements, then the acceleration term may be identified in whole or part with an increment of kinetic energy. For spinning bodies, an increment in the angular velocity will allow the acceleration term to be identified as an increment of the kinetic energy, as is now proved. First, for any arbitrary velocity field, it is possible to consider it as a rigid spin plus an additional part with no angular momentum, as follows. At any time t, the angular momentum of a body about the center of mass is  H = (r × v) dm. (2.23)

9 Some authors use the symbol δ¯ rather than δ to denote this distinction between an increment and the first variation of a functional. 10 (Truesdell and Toupin, 1960, Eq. (218.9)).

Equilibrium figures of spinning solid bodies

The moment of inertia tensor is a functional of the configuration at that time as    J= (x · x)I − x ⊗ x dm. (2.24) From these an “average” angular velocity ω of the body can be defined so that H = Jω

(2.25)

and that defines a spin velocity and acceleration field

as11

vs = ω × x,

(2.26)

as = ω × vs = ω × (ω × x).

The actual velocity and acceleration fields can then be decomposed as v = vs + vˆ ,

(2.27)

a = as + aˆ .

(2.28)

Then it is easy to show that  (x × vs ) dm = H

(2.29)

so that  (x × vˆ ) dm = 0.

(2.30)

This decomposition of the velocity field induces one on the kinetic energy. The kinetic energy is given as  1 KE = (2.31) (vs + vˆ ) · (vs + vˆ ) dm. 2 But since the cross term in the expansion is given as    vˆ · vs dm = vˆ · (ω × x) dm = ω · (x × vˆ ) dm = 0 (2.32) the kinetic energy is just the sum of two positive definite parts:   1 1 ˆ KE = (vs · vs ) dm + (ˆv · vˆ ) dm = KEs + KE. 2 2 (2.33) The spin kinetic energy can be written either as 1 KEs = ω · H 2 or as   1 KEs = (vs · vs ) dm = V dm, 2

(2.34)

(2.35)

where the scalar V is the specific kinetic energy of the spin given as V=

 1 1 vs · vs = ω2 x2 − (ω · x)2 . 2 2

(2.36)

11 The subscript s is used to denote the association with a steady spin.

277

It is easily verified that grad[V] =

∂V = −as ∂x

and  ∂V dm = H. ∂ω

(2.37)

(2.38)

Now, from any initial state, consider any virtual displacements δu; and, for now, also any independent virtual increment δω of the global spin defined above. The increment of spin defines a specific increment of velocity. The variation of the kinetic energy, from Eq. (2.35), using Eqs. (2.37) and (2.38), gives  as · δu dm = −δKEs + H · δω (2.39) but also, using Eq. (2.34)  H · ω. as · δu dm = δKE s − δH

(2.40)

To impose an arbitrary real increment of angular velocity to a real body, it would generally be necessary to apply torques as required from the balance of angular momentum. That is not possible without surface forces. In that case, and for central body forces, the angular momentum H is constant in any actual increment of the motion, and that condition defines a unique permissible actual increment of angular velocity, as can be found by setting the variation of Eq. (2.25) equal to zero: δω = −J−1 δJω.

(2.41)

Then Eq. (2.40) becomes  as · δu dm = δKE s ,

(2.42)

and, assuming conservative body forces and no surface forces, the virtual work Eq. (2.9) gives  δSE + δPE + δKEs + aˆ · δu dm = −δSW d . (2.43) This is the form of virtual work for any increments of configuration and spin related so as to make the angular momentum constant. In the case of an initial state with only an initial spin and no further velocity, the acceleration term is omitted. The application here always assumes that is the case. More generally, without some relation between the general virtual velocity increments and the virtual displacement increments, the remaining acceleration terms are unrelated to a variation of kinetic energy. Finally, along any actual motion, the equivalent of Eq. (2.43) is simply Eq. (2.22) with the kinetic energy decomposition and zero surface traction: d d d [KEs + PE + SE] + KEˆ = − SW d . dt dt dt

(2.44)

278

K.A. Holsapple / Icarus 172 (2004) 272–303

2.3. Average stresses Most often the stress distribution in a body depends not only on its shape and loading, but also on its material behavior. In the study of statics, problems are called “statically indeterminate” when the equilibrium equations and stress boundary conditions do not by themselves suffice to determine the stress state. Then information about the material behavior is needed to derive stress states. However, even though the complete stress state is usually statically indeterminate, certain average stresses in any body are always statically determinate. Of course, the average traction on any coordinate plane cutting across the body is easily found by integrating the body force and acceleration over the cut-away portion of the body, and dividing by the area of the plane. A less well-known result for an average stress for the entire body can be derived, for both statics and dynamics, independent of assumptions about material behavior.12 Starting again with the equation of motion, Eq. (2.7), take the outer product with the position vector x and integrate over the body. Using a Cartesian component index notation,   ∂Tij + ρbi − ρai xk dv = 0. (2.45) ∂xj Vol

Again using the divergence theorem on the first term, one obtains a result for average stresses:    1 1 ¯ Tik dv = ti xk da + (bi xk − ai xk ) dm , Tik = V V i Vol Vol (2.46) where the volume averaged stress tensor T¯ik is the integral of the stress tensor over the body divided by the current volume V . Equation (2.46) then gives the average stresses in terms of the surface, body and inertial forces. Thus, in the static case where the inertial terms are absent, the average stresses in any body are indeed statically determinate; i.e., they are determined from the surface and body loads, without regard to deformation or material properties. In the dynamic case one must know the inertial forces. That generally makes it of less use in dynamics than in statics unless the acceleration field is known; but that is the case for a strained steadily spinning body. This theorem will have important applications below.

3. Stability of equilibrium There are many approaches to definitions and study of the stability of dynamical systems, mostly for special cases. Most are phrased in terms of “infinitesimal stability,” and are 12 For the static case, see Truesdell and Toupin, 1960, Section 220. This is also one of what Chandrasekhar (1969) calls the virial equations of second order.

based on considerations of small perturbations from some initial state. Then there are two further types. In the first case magnitudes of perturbational time histories of a body are considered. In finite-dimensional dynamics, stability is expressed in terms of properties of trajectories in phase space. These cases may require the ability to solve for a time history (or find a Liapunov function). That is also the case when stability is based on the solution for small oscillatory perturbing motions, which is called “spectral stability.” In studies of stability of boundary layers of viscous fluids, the flow is written as a steady mean flow, plus perturbations, then the history of those perturbations are analyzed by the equations of motion. In the second case, some small perturbation at the present time t, most often in the deformation and velocity fields, and some measure of the body with the perturbation at that time t is used directly in the definition, without actually solving for a history. This approach is more common in complex nonlinear theories, but a question arises: can any measure determined only by the present state and certain specific perturbations suffice to determine the subsequent dynamic stability of a system? In the history cases mentioned above, when a linearization is allowed, the stability question often reduces to the investigation of eigenvalues of certain matrices evaluated at the present state. That is, the conditions of the present state do determine the nature of the solutions. That is a consequence of the fact that the governing equations are linear, constant coefficient, second-order differential equations and the increments in deformation and velocity give the required initial conditions to determine the subsequent motions. What in general problems can be determined about boundedness of the motions given only the problem characteristics and “initial” conditions? What “initial conditions” or perturbations need to be considered to determine the subsequent motions? Classical approaches can give some guidance. In the statics of continuum bodies, the “state” is just a geometrical configuration of a body, and small perturbations in configuration only (i.e., displacement) are considered. There cannot be considerations of time or motion. It must be assumed that the configuration determines the stress tensor so there are no rate effects. An example is furnished by the theories of infinitesimal stability of nonlinear elastic bodies (T&N, Section 68bis). They often further assume that the body forces and surface forces are invariant to the increments of displacement.13 For finite-dimensional dynamical systems, stability is based on considerations of increments of both configuration (displacements) and velocities, and the motion resulting from that new “initial state.” In those cases, the “state” consists of a vector of both positions (or generalized coordinates) and velocities (or generalized velocities) of the particles. That state specification determines all forces. An 13 Called “dead loadings,” which includes the absence of loadings.

Equilibrium figures of spinning solid bodies

“equilibrium point” in that state space is a static state only, with zero generalized velocities, but since those can be measured relative to a rigid rotating frame, steady rotating systems are included. Liapunov stability is based on limitations of the extent of the trajectory in the state space from some initial state sufficiently near some equilibrium point. See, for example (Meirovitch, 1997, p. 154). The more restrictive “spectral stability” mentioned above is based on the existence of real frequencies of small oscillations about an equilibrium state. The solutions are found by a linearization of the equations of motion about that equilibrium state. Then, the equations of motion are of the second order, constant coefficient form, so the stability is again determined by the eigenvalues of the system matrices, and is determined by the problem parameters at an equilibrium state. Chandrasekhar (1969) uses this approach for spinning fluid bodies. Here the interest is in a theory of the stability of bodies with self-gravity and a steady-state spin motion, formulated within the context of general continuum mechanics. A general concept of stability would seem to require consideration of all “possible” perturbations from some “state” of equilibrium. The perturbations will cause variations in body forces, inertial forces of spin and internal stresses. Clearly the definition of the “state” of equilibrium needed here is not only a static state, but must include steady motions of a body that maintain a fixed external shape. In Riemann’s ellipsoids, those motions include a steady spin and a steady and constant internal vortical flow. In Maclaurin spheroids and Jacobi ellipsoids, the motion is only a rigid spin. In that case a “state” is a steady strained configuration with a constant spin vector. The definition of the “possible” perturbations also requires discussion. The concept is that of some small but possible variation from an equilibrium state, maybe over some small time-increment. Those increments are more restricted than the virtual displacements and velocities allowed in the virtual work or power statements. That is, the perturbations considered for stability considerations must be consistent with the balance equations, the constitutive equations, and the boundary conditions. I shall call such increments “permissible.” While that requirement is implicit in most studies it is not always explicitly stated nor needed. It is implicit in the methods that use the equations themselves to solve for perturbation histories. It is always assumed that an increment in mass density is related to increments in strain so that the mass elements are constant, i.e., mass balance is enforced. In the static theory it is often necessary to rule out certain rigid rotations from allowable increments, but those can also be ruled out by global balance of moments. In the spectral approach for incompressible fluids by Chandrasekhar (1969) and others, the motion increments are always isochoric; that is a restriction of the incompressible constitutive equation assumption. Then solutions for those increments are found from the balance equations and boundary conditions.

279

So, only perturbations that are consistent with all of the governing equations are to be considered. For an elastic– plastic body at a limit yield state, the cases of “loading” and of “unloading” use different constitutive equations, so that a state of equilibrium at yield is at a state of discontinuous slopes, and that rules out the linearization methods of spectral analysis at those states. Finite-dimensional dynamical systems usually require for their solution initial values only for the displacement and velocity parameters (the state vector in state space). At a static state, perturbations in those two sets, together with the equations of motion, suffice to determine the resulting motion, and consequently, the stability. For the general case, a stability of some equilibrium state at the present time t is to be measured. A displacement perturbation field δu(x, t) over some time δt from the present might be considered. The question is, how arbitrary can this increment be, if it is constrained by the various balance and constitutive equations? It clearly must be consistent with any displacement or velocity boundary conditions. This time dependent perturbation field determines the deformation gradient history and all of its time derivatives, and the velocity field and the acceleration field to time t + δt. It is assumed that the body force field is determined from the geometrical configuration to time t + δt. Balance of mass and the configuration determine the mass density distribution to t + δt. The deformation gradient history over the time interval must be consistent with any constraints from the constitutive equations, and then determines the stress tensor14 and surface tractions. Then, the velocity field and mass distribution determine the angular momentum, which must be consistent with the balance of angular momentum using the body forces and surface tractions. The acceleration must be consistent with balance of linear momentum using the new stresses and body forces. Therefore, the function δu(x, t) cannot be prescribed arbitrarily, but must be constrained by a consistency between boundary conditions, the constitutive equations, and balance of linear and angular momentum. Generally those constraints may be interwoven in a complex way. However, for special constitutive equations, these constraints uncouple and a class of permissible perturbations can easily be delineated. For a purely elastic material, the stress tensor is independent of the history, and depends only on the deformation gradient. Then suppose an arbitrary δu0 (x) is prescribed at time t +δt, subject only to constraints of the constitutive equations and any boundary conditions on displacement. Then in this case that field alone determines the deformation gradient tensor, stress tensor (since elastic), surface tractions and their moment, mass distribution and body forces and their moment at t + δt. Then there are many velocity fields at that time that satisfy the balance of angular momentum, including an instantaneous velocity field of 14 The assumption of a simple material suffices here, but this is true even for “nonlocal” theories.

280

K.A. Holsapple / Icarus 172 (2004) 272–303

a global rigid spin, without affecting any of the quantities already determined. Finally, a consistent acceleration at time t + δt can be found from balance of linear momentum. Therefore, for an elastic material, given any δu0 (x) consistent only with constraints from the constitutive equations and displacement boundary conditions, there is a history δu(x, t), equal to δu0 (x) at t + δt, which is consistent with all equations of the theory. The velocity field at t + δt can be prescribed as any that satisfies balance of angular momentum, for example a rigid spin. The acceleration field can then be solved for using the equations of dynamic equilibrium. Thus, in dynamic elasticity theory, stability considerations can consider any displacement increments, which might just as well be considered instantaneous. For rate-independent materials (e.g., elastic–plastic), it is also true that the stress tensor is determined15 by the field δu0 (x) prescribed at time t + δt, so the displacement increment can be any that satisfies any displacement boundary condition and constitutive equation constraints. For a viscous material with a dependence on the stretching tensor, one can prescribe arbitrary instantaneous displacement and velocity field increments, subject only to the requirement of balance of angular momentum. That can always be satisfied by superimposing a global rigid spin on a prescribed velocity field without affecting the stress tensor. However, if the constitutive equations involve the second or higher-order time derivatives of the deformation gradient, or has more general history dependence, the constitutive equations and the balance of linear momentum again are intertwined, and it is difficult to quantify possible incremental histories in a general way. In that case, while it may be difficult to characterize all possible deformation and velocity perturbation fields, one can still check for instability by considering only a subset of those that are possible. (See below.) Here I will assume that any arbitrary displacement field increment can be prescribed, with a related rigid spin increment satisfying balance of angular momentum. A further velocity field increment is of no consequence for the rateindependent materials assumed here. Stability studies routinely test only a finite subset of all possible perturbations. In the classical work of Riemann, perturbations were limited to constant linear mappings preserving the ellipsoidal figure. In order to study transitions to nonellipsoidal shapes, Poincaré allowed perturbations in ellipsoidal harmonics of higher, but still finite order. Chandrasekhar initially allows arbitrary spatial fields with oscillatory time dependence, but then applies his virial equations of some small finite order, which can only discern low-order polynomial spatial dependences. Fasso and Lewis (2001) allow only perturbations within the “Dirichlet motions” which are also the steady motions determined by a single timedependent second-order tensor (linear mapping determined 15 For sufficiently smooth histories, and small δt, the history can be considered to be of one of the loading, neutral loading or unloading type for the entire interval.

by a 3 × 3 matrix). The velocity field of this mapping is assumed to be steady. Lebovitz (1998) considers with ellipsoidal harmonics through order five. Clearly, if a definition of stability is violated for a set of possible perturbations, it is also violated for any larger enclosing set. Thus, proving a state is unstable using a limited class of perturbations suffices for proving general instability, but proof of stability for a limited class does not prove stability for a larger (e.g., infinite) class. Therefore all conditions of the researchers mentioned above suffice for instability, but are only necessary, not sufficient, for stability. The same will be true here. 3.1. A definition of stability I now consider initial states of an ellipsoidal body with given global spin vector ω superposed on a steady strain field. From that initial “state” I consider infinitesimal increments δu(x) in the displacement field with the unique increment δω in the spin vector that conserves angular momentum, as was given in Eq. (2.41) above. The initial acceleration is only that of the rigid spin, and is denoted by as . The change in the spin acceleration is δas = δω × (ω × x) + ω × (δω × x) + ω × (ω × δu). (3.1) The increment field δu will generate, at each particle, an increment in the body force (self-gravity), an increment in the stress tensor from the constitutive equations, and then Cauchy’s law will determine the accelerations. The acceleration of that perturbed state will not just be a spin, so it can be written as δa = δas + δ aˆ ,

(3.2)

where δ aˆ is the increment of acceleration in addition to that of a global rigid spin with a zero increment of angular momentum. The theorem of virtual work (a consequence of dynamic equilibrium) holds here in the initial state with a = as , and for the perturbed quantities with a = as + δas + δ aˆ . Apply Eq. (2.11) to both states and take the difference to get    δTR : δHT dV − δb · δu dm + δas · δu dm   + δ aˆ · δu dm = δtR · δu dA, (3.3) which can also be written (see Appendix A, Eq. (A.3)) as    δ TR : δHT dV − δ b · δu dm + δ as · δu dm   + δ aˆ · δu dm = δ tR · δu dA. (3.4) Since there is a potential function for the body force and no H = 0, using Eqs. (2.18) and (2.42), surface loading and δH Eq. (3.4) becomes  δ(δSW + δPE + δKEs ) + δ aˆ · δu dm = 0. (3.5)

Equilibrium figures of spinning solid bodies

This general result is now used to propose a measure of stability for a steadily spinning body with steady strain: a state with  δ aˆ · δu dm < 0 (3.6) for all admissible δu(x) (using the steady spin increment defined in Eq. (2.41)) is defined to be stable. From Eq. (3.3) this condition also can be stated as:    R T δT : δH dV − δb · δu dm + δas · δu dm  > δtR · δu dA (3.7) or using Appendix A, Eq. (A.5)  δ

TR : δHT dV − δ  > δ tR · δu dA.



 b · δu dm + δ

as · δu dm (3.8)

For the present case, using t R = 0 and Eq. (2.42), Eq. (3.8) is δ 2 SW + δ 2 PE + δ 2 KEs > 0

(3.9a)

or using the decomposition above into elastic plus dissipative stresses, δ 2 SE + δ 2 PE + δ 2 KEs > −δ 2 SW d

(3.9b)

for all increments of displacement δu(x) consistent with the constitutive equations, balance equations and boundary conditions. Other forms including those with the Cauchy stress are developed in Appendix A. If, in Eqs. (3.9a) and (3.9b), the (=) sign holds for some increments, the state is called neutrally stable, and if the < sign holds for any increment, it is called unstable. Therefore, this is used to define the stability of a spinning equilibrium state with superposed strain. It is a state for which every permissible perturbation in configuration (therefore having the angular momentum balanced spin perturbation), would result in the average increment of acceleration, which is in addition to that of a change in steady spin, being opposed to the increment in displacement. This is a global measure of the tendency to return to a new uniform spin state having the same angular momentum as the unperturbed state. There are many definitions of stability in the literature, and they often can be put into forms involving second variations in energy terms. However, the general forms of Eqs. (3.7) and (3.8) given here apply to a general steady spin state with strains; without any assumptions of statics, or elasticity, or rigidity; or that the body forces are either fixed or conservative, or that the surface forces are fixed or absent.

281

3.2. Stability in special cases Some special applications of this definition of stability are now examined, to show that it gives the expected result. Of course, that does not prove that this approach is ‘correct,’ but does serve to show that it gathers different cases under one umbrella. 3.2.1. Small oscillations about a static state, spectral stability This is the classical free vibrations problem of a deformable body. Suppose the initial state is static with no spin, and assume there is an admissible time-dependent perturbation δu = U(x) Sin(σ t) for some vector field U(x). The present measure of stability considers any δu(x) so we can consider this time-dependent increment at any time t. Notice that for any static state without initial spin, ω = 0 Eq. (2.41) gives δω = 0 so that δas = 0. And, if this is indeed a permissible solution, then the acceleration is at each t, δ aˆ = δa = −σ 2 U(x) Sin(σ t), so directly using Eq. (3.6) gives   δ aˆ · δu dm = −σ 2 Sin2 (σ t) U · U dm (3.10) for any U(x). Therefore, by the present definition, if the body is stable at a static state, the integral on the left is negative so that σ 2 must be positive. Therefore every possible frequency of vibration at a stable state is real. That is what is called spectral stability. Then, stable here implies spectral stability. If there is any solution with imaginary frequency, the state is unstable: unstable spectrally implies unstable here. However, it cannot be asserted that all frequencies being real implies stability according to the present definition, because the test has not been made for all admissible δu(t), but only admissible harmonic solutions. That is, spectral stability does not imply the present stability measure. More importantly, no proof of even the existence of modes of vibration has been proffered, that will depend on the constitution of the material. Generally, in order to satisfy both the equations of dynamic equilibrium only certain U(x) will be possible (eigenvectors). It is only proved here that if there are vibratory motions, according to the present definition of stability they must have real frequencies. 3.2.2. Static strained elastic body In the study of statics in solid mechanics, definitions of stability often restrict material behavior, not body states and motions, although some relation between the two is expected. Consider the dynamic criteria of Eq. (3.7) applied to a static state16 of a strained nonlinear elastic body, and first assume that the body forces and surface forces are unchanged by displacement increments.17 All spin terms can 16 In a static theory, only the constitutive equations for static strains need be specified. Here the assumption is that the material is elastic for all strain histories. 17 The “dead loading” mentioned above.

282

K.A. Holsapple / Icarus 172 (2004) 272–303

be dropped from all expressions. In this case, Eq. (3.7) reduces to the single term  δTR : δHT dV > 0 (3.11) the form given by T&N, Eq. (68b.11) as their definition of the “superstability” of a general elastic body in statics. If the > sign is replaced by , they call the state simply “stable” as in their Eq. (68.b1), here that includes neutral stability. Since Eq. (3.7) applied to a static state implies a dynamic acceleration increment satisfying Eq. (3.6), it is proved that the material requirement of Eq. (3.11) is indeed the same as the dynamic consideration of Eq. (3.6). It should be noted that many other material criteria have also been studied, with a similar relation but using different stress and deformation measures. Those represent different physical requirements (see Lubarda, 2002, p. 165). T&N go on to assume that the increment in stress can be given in component form as δTR ij =

∂  R T δHmn = Aij mn δHmn ∂Fmn ij

(3.12)

using the “instantaneous elasticities” (in the strained state x) as Aij mn =

∂TR ij ∂Fmn

.

Then the Eq. (3.11) becomes (T&N, Eq. (68b.11))  Aij mn δHj i δHnm dV > 0

(3.13)

(3.14)

Vol

for all δH. Hadamard’s theorem replaces this global statement with a local requirement on the elasticities at each point in the body (see T&N, Eq. (68b.3)). The essential idea of that proof is that the integral inequality must hold for all fields δH, and one that is nonzero only locally around any point can be used. Further, T&N prove that the local condition is sufficient to ensure uniqueness of solutions. And, they prove that with zero surface forces, every possible frequency of infinitesimal free sinusoidal motion about a stable configuration of an elastic body is real and when the stronger (>) sign holds every frequency is positive. Thus, for an elastic body with dead body loads, the stability definition here again implies spectral stability. The converse, that spectral stability implies (general) stability is not true. T&N consider also for a purely static case the more special case with a potential energy and a hyperelastic material. They define the state function J of the configuration as, using the notation here, (T&N, Eq. (88.10))    J(X) = Σ(F) dm + U (x) dm − u · tR dA. (3.15) Further, they assume that everywhere on the boundary either δu = 0, or δtR = 0 so that the last term has zero variation. Then, their definition of stability is that δ 2 J  0, which, is just Eq. (3.9a) here without the kinetic energy term, and with

the stress work as a strain energy increment. If δu is not zero on the surface and the potential term is absent, the surface term remains, and that expression leads to stability of static elastic systems in structural mechanics such as classical (Euler) beam buckling and other structural buckling problems. Note that an underlying assumption made above is that the solutions about an equilibrium state are smooth, so the linearization used in Eq. (3.13) is valid. The linearized problem admits oscillating solutions. Again I note that plasticity solutions do not have that smoothness, which rules out using these approaches. In many cases an increment in one direction is not even admissible in the opposite direction, as will be discussed below. 3.2.3. Static elastic–plastic materials Next, for an elastic–plastic material in a static state, the inertia and acceleration terms are also absent. Again assuming “dead ” body force terms,18 Lubarda (2002), (Eq. (10.6.4)) also defines the stability of a body of a general elastic–plastic material with the equivalent19 of Eq. (3.11), but includes a term with an increment of surface traction as in Eq. (3.7). He attributes the definition to Hill (1958, 1978). He proves that the incremental solution need not be unique at a stable state, but that bifurcations can occur as long as the body and surface forces are not dead loadings. In plasticity theories, a common definition of material stability is called Drucker’s postulate. Chen and Han (1988), (Eq. (3.158)) give that postulate, for the infinitesimal strain and static case, as the condition that   δb · δu dm + δt · δu da > 0 (3.16) for any part of a body; and further that such an expression is satisfied also for cycle of addition and removal of force increments. They ignore the difference between deformed and undeformed volume and surface integrals. Using this in the static case of Eq. (3.7) gives that the stress integral must be positive, but now for every volume. Therefore one gets the material condition20 δT : δe > 0

(3.17)

for every strain increment δe. In one-dimensional stress test this implies that the stress-strain curve always has positive slope, a plasticity material property called strain hardening. Lubarda (2002) develops the theories of plasticity in a complete nonlinear, large-deformation framework. He states an equivalent of Drucker’s postulate in that nonlinear framework, and proves that it implies an invariant equivalent to 18 It is obvious that if he admitted nonzero body force increments, Lubarda would also include the body force term. 19 His “nominal” stress tensor is the transpose of the Piola tensor used here. 20 Chen and Han (1988) state this in terms of time derivatives of a possible motion, not increments. That is equivalent to increments that satisfy the balance equations. The context is the infinitesimal theory only.

Equilibrium figures of spinning solid bodies

Eq. (3.17) that uses the Piola stress and a material plastic strain tensor, his Eq. (8.6.10). He also discusses another postulate called Ilyushin’s postulate, based on a dual formulation where the stress tensor is decomposed into an elastic and plastic parts. These are clearly stronger requirements than simply requiring the sum of terms of Eq. (3.7) integrated over the entire body to be positive. It is a special assumption about material behavior, and plastic materials that do not satisfy these conditions are also considered. The use of Drucker’s postulate leads to “associated” flow rules, otherwise the flow rules are called “nonassociated.” Below features of plasticity called “strain softening” are considered, and those cases do not satisfy Eq. (3.17). 3.2.4. Liapunov stability In the theory of the dynamics of systems with a finite number of degrees of freedom, the definition of Liapunov stability is well known and studied. Here the formulation has an infinite number of degrees of freedom. However, as mentioned above, stability tests are routinely restricted to a finitedimensional subset of all possible increments. Therefore, an effective reduction to a finite-dimensional case occurs. For some finite-dimensional set of virtual displacements that conserve angular momentum, referring to Eq. (2.44) above, we can define the dynamic potential as U = SE + ˆ so that PE + KEs , and the Hamiltonian21 as H = U + KE, Eq. (2.44) is d d H = − SW d . (3.18) dt dt If Eq. (3.9b) is satisfied and the dissipative stress work is positive semi-definite and never decreasing, then U is posˆ is positive definite (see Eq. itive definite; and, since KE (2.33)), H is a positive definite function, with a nonpositive time derivative. Thus, a state of a general body that is stable by the present definition with nonnegative dissipative stress power will be stable for finite-dimensional sets of virtual increments according to the Liapunov criteria. (See Meirovitch, 1997, p. 157.) Fasso and Lewis (2001) use such a Hamiltonian in their study of the stability of the Riemann ellipsoids for the fluid case, where the dissipative stress work is exactly zero.

4. Soil behavior, plastic yield conditions, and flow rules The present application is to bodies presumably composed of the same “geological” materials as on earth, especially soils and rocks. Those are known to be highly nonlinear. The elastic response from any given state depends significantly on that state. For example, data for dense dry sands show that the shear modulus increases by a factor of four 21 Compare to Meirovitch (1997) (Eq. (2.160)), for the finite-dimensional equivalent.

283

as the confining pressure changes from 24 kPa (the average gravitational pressure in a spherical body of density 2 g/cm3 and radius 10 km) to 216 kPa (the average in a 30 km body). That modulus also decreases by a factor of almost 2 as the void ratio increases from 0.35 to 0.75. Repeated loadings can greatly diminish the shear modulus. Even when elastic, these materials are certainly not linear elastic. An assumption of linear elastic behavior is not appropriate and will not generally be used here, except as a local approximation for some small increments from an arbitrarily strained state. Similarly, the plastic yield behavior can be complex. In tests with increasing strain, the stress after a peak can decrease markedly. See for example (Lambe and Whitman, 1969, Chapter 10). This phenomenon is called strain softening, and is a common feature of geological materials. Even more dramatic strain softening is apparent in the curves for concrete given in (Chen and Han, 1988, Chapter 7). This feature of strain softening leads to possibilities of strain localization and other material instabilities. 4.1. Perfectly plastic theory For the description of the initial yield limits, a common assumption for soils and rocks is that the limit states are governed by a Mohr–Coulomb yield criteria22 (see, e.g., Chen and Han, 1988). While any of the stress tensors may be used to define yield functions (see Lubarda, 2002), standard plasticity models use the Cauchy stress tensor and use the symbol σ for that Cauchy stress. I shall now follow that convention. The Mohr–Coulomb model is based on the assumption that there is a limit shear stress magnitude τ on any plane in the body that increases linearly with the compressive normal stress σn (a negative quantity) on the plane: τ = Y − σn Tan(φ),

(4.1)

where the coefficient of proportionality is the tangent of the angle of friction φ. This defines an envelope of maximum shear stress, and the Mohr’s circles of all possible stress states must lie within that envelope. Figure 1 shows the yield envelope with a typical Mohr’s circle: It is convenient to use two different subscript labeling schemes. The x, y, z subscripts will denote the components in those coordinate directions. The subscripts 1, 2, 3 are used for the directions of the principle stresses. If σ1 is the 22 Students of soil and rock mechanics will know that there are two different but similar yield functions used to describe pressure-dependent yield. The Drucker–Prager relation is algebraically more complex, but a single algebraic form holds for all stress combinations. For presentation of definitive and simple formulas such as that for the maximum spin in Eq. (6.3) below, the Mohr–Coulomb form is more convenient, but six different combinations of the ordering of principle stresses occur and must be considered. For use in a symbolic program such as “Mathematica” the Drucker–Prager form is more useful but gives more complicated formulas. Both approaches give very similar results. Here I have chosen to present the analysis and results using the Mohr–Coulomb form, although I have done it both ways.

284

K.A. Holsapple / Icarus 172 (2004) 272–303

Fig. 1. The Mohr–Coulomb yield criteria. The quantity Y , the maximum shear stress at zero normal stress, is called the cohesion. The maximum normal (tensile) stress is σtensile = Y /Tan(φ). When the cohesion Y is zero, tensile stress is not allowed, but shear and compressive strength remains.

largest and σ3 is the smallest principal stress, then the Mohr– Coulomb yield criteria is given for general stress states by

(σ1 − σ3 ) 1 + Tan2 (φ) + Tan(φ)(σ1 + σ3 )  2Y. (4.2) This can also be written as Y mσ1 − σ3  2 ( 1 + Tan (φ) − Tan(φ))

(4.3)

with m=

1 + Sin(φ) . 1 − Sin(φ)

(4.4)

On a three-dimensional plot of the three principle stresses σ1 σ2 σ3  this defines a limiting cone with a regular hexagonal cross-section which increases in size as the negative of the sum of the three stresses (the pressure) increases. The intersection of this surface with a plane of constant pressure σ1 + σ2 + σ3 = constant is depicted in Fig. 2. For rubble-pile materials, the cohesion (and therefore the tensile strength) is assumed to be zero. Then the apex of the hexagonal cone is at the origin in the above figures. The yield criteria is then given simply as23 mσmax − σmin  0.

(4.5)

In this case the material strength is determined solely by the angle of friction φ (or the parameter m from Eq. (4.4). To completely describe a plastic model, one must also prescribe the “flow rule” which governs plastic strain increments when yielding (loading), and the behavior for nonyielding (unloading) increments. The plastic flow rule can either be “associated” or “nonassociated.” The associated rule has the plastic strain increment during plastic flow (loading) normal to the yield function surface. 23 Note that the stresses are in compression, and therefore are negative. The smallest has the largest magnitude, and is the most compressive.

Fig. 2. The yield surface for a Mohr–Coulomb material in a plane of constant pressure σ1 + σ2 + σ3 = constant. Strain increments are shown at a point where the largest principle stress is σ1 and the smallest is σ3 , so the surface is given by mσ1 − σ3 = constant. For a state on the yield surface and “loading,” the plastic strain increment is defined by the flow rule and always in the outward direction as shown. For the associated flow rule it is normal to the yield surface. For “unloading,” any of the increments pointing inward from the surface are possible, and will be elastic. Other directions of strain increments are not admissible.

That is in fact a consequence of Drucker’s stability postulate mentioned previously. In that case, for any plastic yield function of the type f (σ1 , σ2 , σ3 ) = 0, the plastic strain increment is defined by the form δep = λ

∂f gradσ f, ∂σ

(4.6)

where  > 0 if f = 0 and grad f · dσ > 0 (loading), σ    = 0 if f = 0 and grad f · dσ < 0 (unloading), σ λ=  = 0 if f = 0 and grad  σ f · dσ = 0  (4.7) (neutral loading). These are for a work-hardening material. For a perfectly plastic material the first case is absent, and the plastic strain is nonzero in the last case. For the nonassociated case, a form similar to Eq. (4.6) is commonly assumed, but using a potential function g rather than the yield function f . Both cases will be included below in the applications.

Equilibrium figures of spinning solid bodies

An equation for the nonplastic increments is also required, it is commonly assumed to be a linear elastic relation between the increments of strain and stress. Since it is only for the increments, the corresponding elastic modulii can be allowed to be dependent on the stress, thereby including nonlinear elastic behavior. Therefore, in a perfectly plastic body, a strain increment from some current state can be either elastic or plastic. States within the yield surface have only (incrementally) elastic strain increments corresponding to increments of stress. If the stress state is on the yield surface and the stress increment is “unloading,” then the plastic increment of strain is zero and the entire strain increment is (incrementally) elastic. Finally, if on the yield surface and “loading” then the plastic strain increment is nonzero and, from Eq. (4.6), geometrically is perpendicular to the yield surface, as is shown in Fig. 2. For the Mohr–Coulomb criteria, if the components are ordered by directions of the maximum, intermediate, and least principal stress, then the associated flow rule obtained from Eqs. (4.5) and (4.6) can be written in the simple 3 × 1 matrix form (Chen and Han, 1988, Section 4.50) δε p = λm 0 −1

(4.8)

for some positive scalar λ. This is the ordering of the components in the principle stress numbering scheme. The ordering in x, y, and z components may differ. When, for example the largest principle stress is σx and the smallest is σy , the x, y, and z components of strain increment corresponding to Eq. (4.8) are given as λm −1 0. The ordering to use will be apparent from the use. During any plastic flow, assuming again the ordering from largest to smallest principle stress, since the yield function is given as mσ1 − σ3 = 0, the three principle stress components have the form24 σ = σ1 σ2 mσ1 .

(4.9)

The plastic work increment involves the integral of the dot product of the vectors of principal stress and principal plastic strain increment:   m σ : δep = λσ1 σ2 mσ1   0  = 0, (4.10) −1 which is zero for any plastic loading. Thus, the plastic work increment is always zero for a cohesionless perfectly plastic Mohr–Coulomb material with the associated flow rule, since the stress vector to the yield surface is always normal to the plastic strain increment. That is obvious from the above Fig. 1 when the apex of the yield cone is at the origin: then the stress vector is along the cone, and the plastic strain increment is normal to it. This property would not be true if a 24 Here the symbol σ stands for the vector of principle stresses, elsewhere it is a tensor, the context should make clear which it is.

285

nonassociated flow rule was adopted, or for a material with cohesion. 4.2. Work hardening and softening, general flow rules In fact, in this study more general plasticity models will be included. Without giving detailed equations, I shall use the theory given in (Chen and Han, 1988, Section 5.7). The yield function is the Mohr–Coulomb model, but adds possible work-hardening or softening. The work-hardening modulus (Chen and Han, 1988, Eq. (5.140)) κ is a function of the plastic work (or equivalent plastic strain) and takes on positive, zero, or negative values for a material that is, respectively, work hardening, perfectly plastic, or work-softening. The flow rule may be associated or not, but is governed by a loading function. The associated flow rule is perpendicular to the yield function, as given in Eq. (4.8), otherwise it is determined by the same (Mohr–Coulomb) form, but using a different scalar loading parameter m1 with 1 < m1 < m. The case with m1 = 1 is incompressible, the case with m1 = m has bulking (a plastic volume increase when yielding). The incremental elastic strains are for an elastic material, although there need not be a global elastic model. They are defined using at least two independent elastic modulii, a bulk modulus K and a shear modulus µ. Then, in this general case there is, at a given stress state and with a given plastic work, a definite but complicated relation giving the stress increments corresponding to any strain increments, given by (Chen and Han, 1988, Eq. (5.144)). It has the form δσ = Cep δe,

(4.11)

where Cep is the fourth-order elastic–plastic incremental stiffness tensor given by a rather messy expression in their Eqs. (6.104), (5.135), (5.145), and (5.146), which involve derivatives of a yield function and of a loading function with respect to the stress components. In all cases here the strain increments will be for the principal strains only, denoted by the vector δε, and then the tensor expression Eq. (4.11) can be written more simply as a vector equation using a secondorder modulus tensor in component form as ep

δσi = Cij δεj .

(4.12)

The modulus tensor is easily evaluated using “Mathematica” from the equations of Chen and Han, and is determined (in the isotropic case) by 5 material parameters: the elastic modulii K and µ, the yield function parameter m, the flow rule parameter m1, and the hardening modulus function (of the plastic work) κ. The analyses here will use this general plastic incremental relation. 4.3. Material properties of a typical sand To make estimates from this general theory and to understand the meaning of the parameters, one must know how the model relates to experimental tests on actual soils or rocks. For any stress test the theory gives definite results,

286

K.A. Holsapple / Icarus 172 (2004) 272–303

Fig. 3. Triaxial compression of a dense sand. The lateral stresses are at 100 kPa pressures while the vertical stress and strain increase in compression. The results are approximately linear elastic to about 100 kPa.

which can be evaluated using “Mathematica.” Consider a typical triaxial compression test of a soil; in that case compressive lateral stresses σh are held constant and the vertical stresses σv increase in compression. Some experimental results for (σv − σh )/2 versus the vertical strain εv are as shown in (Lambe and Whitman, 1969, Fig. 10.18) for loose and dense sands and in Fig. 10.13 for a dense sand. The experiment of the latter, from Fig. 10.13, had an initial pressure of 100 kPa. The measurements of the shear stress = (σv − σh )/2 are reproduced here as Fig. 3. To fit the model to this data, I assume a Poisson’s ratio of 0.25, an angle of friction of 30◦ so that m = 3, and I take m1 = m. With these values initial yield is predicted when σv = mσh = −300 kPa or τ = 100 kPa, consistent with the data. For this Poisson’s ratio, the model gives the elastic relation as 5 dτ = µ dεv , (4.13) 4 where µ is the elastic shear modulus, and the data to the shear stress of 100 kPa gives µ = 30 MPa. Above τ = 100 kPa, the model gives for the plastic increments   5r µ dεv , dτ = (4.14) 10 + 4r where r = κ/µ is the ratio of the work-hardening modulus to the shear modulus. The work-hardening modulus function κ can be found to exactly fit the data points, and the result is shown in Fig. 4. The Fig. 10.18 in (Lambe and Whitman, 1969) shows a much larger negative slope and work softening for a different dense sand. The conclusion is that the work-hardening parameter can easily become a negative value, indicating work

softening, and in magnitude is equal to a few percent of the shear modulus value. 4.4. Large strain theory One final topic must be addressed. The plasticity theory of Chen and Han is only for the infinitesimal theory. It uses the Cauchy stress tensor and the Eulerian strain. It is well known that a proper theory for large strains with general functional dependences cannot use the Cauchy stress and Eulerian strain. Lubarda (2002) presents the nonlinear theory of plasticity. He shows that an incremental elastic–plastic stress-strain law holds for various higher-order theories using what he terms “conjugate” stress and strain measures. It is convenient to use the present state as a reference for the measurement of increments. Then he gives the specific form for a stress rate ep τ˚ = L(0) and a strain rate D as (using his notation) in his Eq. (9.4.19) (and in other places) as ep

τ˚ = L(0) : D,

(4.15) ep

where the fourth-order elastic–plastic modulus tensor Λ(0) has an elastic and a plastic part. Using his Eqs. (2.5.13) and (3.9.12) gives the Cauchy stress rate as ep

σ˙ = L(0) : D + Wσ − σ W − (Tr D)σ .

(4.16)

The time derivative of the stress plus the two terms with the spin tensor can be recognized as the Jaumann stress rate. ˙ where With the present as reference, then D = e˙ and W = R, R is the rotation tensor. Therefore the incremental form of Eq. (4.16) is ep

δσ = L(0) : δe + δRσ − σ δR − (Tr δe)σ .

(4.17)

Equilibrium figures of spinning solid bodies

287

Fig. 4. The work-hardening (softening) modulus as a function of the plastic strain from the data of Fig. 3. It begins as hardening with about 40 MPa (off the top of this scale) and becomes softening (negative values) at a plastic strain of about 3%. The maximum softening is a little less than 1 MPa at a strain of about 5–6%. That is about 3% of the shear modulus.

An experiment with no rotation can then be used to determine the elastic–plastic tensor in this relation. A specific form using one set of measures does not translate into the same form using another. Differences are on the order of the stresses compared to the modulus terms, which are generally very small numbers. Here, I shall simply assume that the theory from Chen and Han gives directly the first term of the above relation: that is, I assume that the modep ulus L(0) in Lubarda’s notation is the same as the tensor Cep in Eq. (4.11) above. Thus, the constitutive equation to be used is δσ = Cep : δe + δRσ − σ δR − (Tr δe)σ .

(4.18)

In the experiments presented above, the stresses are several orders of magnitude less than the modulii values. Therefore one could equally well assume Cep included both the first and last terms in Eq. (4.16). However, if any results depend in any essential way on such subtle second-order differences in the model, they would be suspect anyway, since any model is at best an approximation to the actual behavior of a Solar System body. So, using this particular model, I look for general features of the solution determined by the plasticity and plastic deformation.

a b c = ax ay az , the gravitational potential is given (e.g., Chandrasekhar, 1969) in the rigid coordinate system rotating with the body as   U = πρG −A0 + Ax x 2 + Ay y 2 + Az z2 , (5.1) where the A0 and Ai (i = x, y, or z) terms are given as   du du , Ai = abc A0 = abc , (5.2) ∆ ∆(ai2 + u) and ∆ = (a 2 + u)(b 2 + u)(c2 + u). They are related by Ax + Ay + Az = 2, 2

2

(5.4)

Denoting the ratios of the semi-axes lengths as c b α= , (5.5) β= a a the three Ai can be put into the nondimensional forms: ∞ Ax = αβ 0

∞ Ay = αβ 0

5. Spinning ellipsoidal bodies The theory is now applied to the application of interest, spinning solid and fluid ellipsoidal bodies with self-gravity.

(5.3)

Ax a + Ay b + Az c = A0 . 2

du , (u + 1)3/2(u + β 2 )1/2(u + α 2 )1/2 du , (u + 1)1/2(u + β 2 )3/2(u + α 2 )1/2

∞ Az = αβ 0

du (u + 1)1/2 (u + β 2 )1/2 (u + α 2 )3/2

(5.6)

5.1. Gravitational forces

showing that they are dimensionless and depend only on the shape, not the size of the body. It is easily shown that for a spherical body of radius R, that

For a homogeneous ellipsoidal body, with the lengths of the semi-axes labeled in order from largest to smallest as

2 Ax = Ay = Az = , 3

A0 = 2R 2 ,

(5.7)

288

K.A. Holsapple / Icarus 172 (2004) 272–303

but for a general ellipsoidal body the results are not so simple and generally require numerical evaluation. In any case, the body force components are given as ∂U = −2πρGAi xi (with no sum on i). (5.8) ∂xi Using the integrals  2  ai xi xk dv = 5 V , i = k, (5.9) 0, i = k, the potential energy of the body, using Eqs. (5.4) and (2.14) is given as

bi = −

4 PE = − πρGMA0, (5.10) 5 where M is the mass of the body and where V denotes the present volume. 5.2. Kinematics of spin The body is assumed to have a constant spin of rate ω about its shortest principle z-axis, the spin state with the least kinetic energy. The velocity x, y, and z components are v = ω−y x 0,

(5.11)

and the acceleration components are a = −ω2 x y 0.

(5.12)

The component about the z-axis of the angular momentum of Eq. (2.25) is written in scalar form as H = Jz ω,

(5.13)

where the rotary moment of inertia component Jz is about the shortest axis and given for an ellipsoidal body with mass M as  M 2 a + b2 . Jz = (5.14) 5 The kinetic energy of spin is simply 1 KEs = Jz ω2 . 2

(5.15)

5.3. Average stresses in a spinning body First, the area average across any coordinate plane can easily be found by a single integration. For example, for the plane x = x0 , the average normal stress σx is calculated by integrating the x-components of the acceleration given by Eq. (5.12) minus the body force given by Eq. (5.8) over the portion of the body x > x0 and dividing by the area of that plane. The result is given as   2 3/2  a2  x0 σˆ x = ρω2 − 2πρ 2GAx (5.16) 1− 4 a with similar results for the other coordinate planes (the z planes do not have the spin term). The maximum area average over all such cross-sections is on the plane cutting across the center x0 = 0 and is

  a2 σˆ x = ρω2 − 2πρ 2 GAx (5.17) 4 and again the obvious changes give the results for the other two coordinate planes through the center of the body. The more important averages for the present are the volume averages for the entire body. Those can be determined by using Eqs. (5.8) and (5.12) in Eq. (2.46). The surface term is zero since there are zero tractions. Then Eq. (2.46) gives the volume average stresses over any spinning ellipsoidal body as  a2  , σ¯ x = ρω2 − 2πρ 2 GAx 5   b2 , σ¯ y = ρω2 − 2πρ 2 GAy 5  c2  σ¯ z = −2πρ 2GAz (5.18) , 5 and the volume average shear stresses are zero. The relative shape of the body determines the three Ai values from Eq. (5.6). These are the volume average stresses in any spinning ellipsoidal body, whether the body is comprised of rock, rubber or regolith. Every analysis of the stresses or stress resultants using any particular material model (rigid, linear or nonlinear elastic, fluid, viscous, viscoelastic, elastic–plastic, perfect gas, etc.) must have these same volume averages. The negative of the sum of the three, divided by 3, and using Eq. (5.4) gives the volume average pressure in any spinning ellipsoidal body with self-gravity as (a 2 + b2 ) 2π 2 + ρ GA0 . (5.19) 15 15 Note that the area averages across the central coordinate planes are just 5/4 of the volume averages. The actual maximum stresses, at the center, depend on the material model and are not statically determinate. However, Holsapple (2001) shows that the yield condition requires all normal stresses to be zero at the surface of the ellipsoid. It is not hard to surmise that the normal stress distribution across any coordinate plane is quadratic. If it is zero at the surface, using as an example the coordinate plane x = 0, the normal stress distribution on that plane must have the form   2  2 y z − σˆ x = K 1 − (5.20) b c

p¯ = −ρω2

and by equating its average over that plane to Eq. (5.17) it is found that the value of the constant K and the maximum σx is just two times the area average of Eq. (5.17) and 5/2 the volume average of Eq. (5.18). That is just the result in (Holsapple, 2001) where the entire stress distribution was derived. Other solutions such as elastic ones do not have such simple results. 5.4. Virtual work for spinning bodies I now derive the special form of the virtual work and the stability measure for this ellipsoidal body. Following the lead

Equilibrium figures of spinning solid bodies

of others mentioned above, I investigate only displacement increments consisting of a linear mapping of the spatial coordinates defined by a constant tensor δB: δu = δBx.

(5.21)

δ 2 u = δBδx = δB(δBx).

δH = δB,

δ 2 H = δBδB.

(5.23)

The symmetric part of the constant tensor δB defines constant stretches and shears, and the skew symmetric part a rigid rotation increment. These strain increments might be elastic or plastic, or both. For consideration of stability, elastic states admit arbitrary values for this tensor; but plastic increments cannot be arbitrary, and will be required to be consistent with the flow rule. The spin acceleration given in Eq. (5.12) has only x and y components. Therefore the acceleration terms from virtual work are    2 x(δB11 x + δB12 y + δB13 z) as · δu dm = −ω Mass

mass

 + y(δB21 x + δB22 y + δB23 z) dm. (5.24)

The symmetry of the ellipsoidal body makes the cross terms involving products such as xy integrate to zero, so that using Eq. (5.9) for the remaining integrals gives   Mω2  δB11 a 2 + δB22 b2 . as · δu dm = − (5.25) 5 The increment in the stress work is simply     R T σ dv : δe T : δH dV = σ : δe dv = = V (σ¯ x δB11 + σ¯ y δB22 + σ¯ z δB33 )

(5.26)

expressible in terms of the average stresses and the three normal strain increments and the present volume V (recall that the average shear stresses are zero). No constitutive relation is needed to calculate this stress work, only the average stresses. Shear strains and spins play no role. Finally, the body force components from Eq. (5.8) are b = −2πρGAx x Ay y Az z

5H , M(a 2 + b2 )

(5.29)

(5.22)

For the calculations of various increments, it is most convenient to take the present state as the reference. Then, using Eqs. (5.21) and (5.22) other kinematical relations are, to appropriate orders F = I + δB,

Finally, on the surface δt = 0. Because we will be interested in variations with fixed mass M and angular momentum component H , the mass density ρ and angular spin ω are written in terms of those two quantities as: ω=

This has a second-order variation

289

(5.27)

and again, using the fact that integrations of product terms are zero, the integrated body force terms are  Ma 2 b · δu dm = −2πρG 5   × δB11 Ax + β 2 δB22 Ay + α 2 δB33 Az . (5.28)

3M . (5.30) 4πabc Equations (5.25) and (5.28) involve only the normal strain increments, which are now denoted more concisely in a matrix notation as      δB11   δεx  δε = δB22 = δεy . (5.31)     δB22 δεz

ρ=

The change in the lengths of the axes are then given by the vector δs = δa δb δcT = δεx a δεy b δεz cT .

(5.32)

Then, by grouping the various terms the virtual work theorem of the form of Eq. (2.43) can be combined into the single vector form25 δKEs + δPE + δSW = E · δε = 0

(5.33)

with the matrix (vector) E given as T  −5H 2a 2 4πabc  3M 2 GAx a    + σ¯ x +    M(a 2 + b2 )2 10bc 3          3M 2 GAy b 4πabc −5H 2 b2 . (5.34) E= + σ ¯ + y  M(a 2 + b2 )2 10ac 3            3M 2 GAz c 4πabc     + σ¯ z 10ab 3 This vector includes the variation of the kinetic energy of spin (the terms with H ), the variation in potential energy (the terms with G), and the stress work increment (the terms with the average stresses). It might be noted that it has the units of energy, or stress times volume. The spin terms are characteristic of spin stresses, and the gravity terms of gravity stresses, both times the volume of the body. The theorem of virtual work states that at an equilibrium state, Eq. (5.33) must be zero for all increments δε. Therefore, if the average stresses of Eq. (5.18) (which were derived assuming equilibrium with a constant spin) are inserted into Eq. (5.34) one simply finds that the vector E is identically zero; this is nothing more than a check on the algebra. (In fact, this could be used to derive the expressions Eq. (5.18) for the average stresses.) However, for testing for stability, the expression applicable to incremental states 25 Since the matrix E is defined as a row matrix and δε as a column, no sign is needed in this matrix expression. I include the dot-product symbol (·) for those that prefer to think of these matrices as vectors. Also, δSE and δSW d are regrouped into δSW.

290

K.A. Holsapple / Icarus 172 (2004) 272–303



that are not in dynamic equilibrium with a constant spin are needed in order to take the appropriate derivatives. Therefore, it is not appropriate to evaluate the average stresses in Eq. (5.34) before considering the second variations and resulting derivatives. 5.5. Stability of spinning bodies For stability considerations, the second-order variations are needed. For brevity, I shall ignore shear increments and restrict the increments to those of normal strain only. Eq. (3.9a) gives δ(δKEs + δPE + δSW) = δ(E · δε) > 0.

(5.35)

Therefore, the stability criteria requires that δ(E · δε) > 0

(5.36)

for every δε consistent with fixed mass, fixed angular momentum, and all restrictions from the constitutive equations. This vector E will be considered in two parts. The kinetic energy and potential energy terms are denoted as E1 and given as  T −5H 2a 2 3M 2 GAx a     +     10bc  M(a 2 + b2 )2        2 2 2 3M GAy b −5H b E1 = (5.37) . +   M(a 2 + b2 )2 10ac           3M 2 GAz c     10ab This part is determined entirely by the configuration, i.e., by the axes lengths a, b, and c, for given mass M and angular momentum H . The variation of that part is due solely to the variation δs of the axes lengths from Eq. (5.32). Reverting to the index notation ai for these three axes, the contribution of this part to the stability criteria Eq. (5.36) is given as δ(E1 · δε) =

∂E1i δaj δεi . ∂aj

(5.38)

Equation (5.32) gives the increments of axes lengths. If one defines a new tensor G1 as G1ij =

∂E1i aj ∂aj

(no sum on j )

(5.39)

then Eq. (5.38) becomes δ(E1 · δε) = G1ij δεj δεi = δεT G1 δε.

(5.40)

The derivatives of the E1 with respect to the ai needed in Eq. (5.39) can be evaluated from their definitions in Eq. (5.37). (Note that M and H are fixed.) In the process, derivatives of the three potential terms Ai with respect to the axes lengths are needed. Their definitions in Eq. (5.6) can be used to find those derivatives in terms of integrals defined as  du , Ia = 2 3/2 2 (a + u) (b + u)1/2 (c2 + u)1/2

du , Iaa = (a 2 + u)5/2(b2 + u)1/2 (c2 + u)1/2  du Ib = , (a 2 + u)1/2 (b2 + u)3/2 (c2 + u)1/2  du , Iab = 2 3/2 2 (a + u) (b + u)3/2 (c2 + u)1/2  du Ic = , 2 1/2 2 (a + u) (b + u)1/2 (c2 + u)3/2  du Iac = , (a 2 + u)3/2(b2 + u)1/2 (c2 + u)3/2

(5.41)

and similarly for Ibb , Icc , and Ibc .26 With these definitions, the needed derivatives of the potential terms become  bcIa − 3a 2 bcIaa acIa − ab2 cIab abIa − abc2 Iac   ∂Ai = bcIb − a 2 bcIab acIb − 3ab2 cIbb abIb − abc2 Ibc ∂aj bcIc − a 2 bcIac acIc − ab2 cIbc abIc − 3abc2 Icc (5.42) which can then be used in Eq. (5.39) to determine the G1 matrix. 5.6. The strain energy terms To the spin and gravity terms must be added the stress work terms, which depend on the material model, and given  generally as δ TR : δHT dV or δ σ : δe dv. Note however, by the previous remarks, one cannot just use the stresses assuming equilibrium in this relation. Instead, one must use the constitutive equations to calculate the stress increments from given strain increments. The first two terms of Appendix A Eq. (A.11) give the stress terms, using now the symbol σ for the Cauchy stress, as   2 δ SW = δσ : δB dv + Tr[δB]σ : δB dv. (5.43) The general relation between strain increments and stress increments for a general work-hardening elastic–plastic body was assumed as in Eq. (4.18) with four terms. At present, the increment of displacement has a strain increment δB with zero rotation increment δR, leaving two terms. Using those, the increment of the stress work Eq. (5.43) is derived very simply27 as   2 δ SW = δσ : δB dv + Tr[δB]σ : δB dv   ep  = C δB − Tr[δB]σ : δB dv  + Tr[δB]σ : δB dv

26 These are essentially the same as those denoted as A by Chandraij sekhar (1969). 27 The astute reader will recall the discussion following Eq. (4.18), and note some arbitrariness regarding the Tr[δB] terms.

Equilibrium figures of spinning solid bodies



 =

(C δB) : δB dv = ep

 C dv δB : δB ep

¯ ep δB : δB), = V (C

Define an average radius as (5.44)

where, recognizing that the elastic–plastic modulus terms are not generally constant throughout the body, their vol¯ ep . Then, using ume average is indicated by the symbol C the normal strain increments of Eq. (5.31), the last term in Eq. (5.44) can be written using the elastic–plastic modulus 3 × 3 matrix of Eq. (4.12) and the three normal strain increments as ep δ 2 SW = V C¯ pq δεp δεq .

(5.45)

Defining the symmetric matrix G2 as 4π ¯ ep abcC ij 3 then Eq. (5.45) is just (G2 )ij =

δ SW = δε G2 δε 2

T

291

(5.46)

r¯ = (abc)1/3

(5.51)

and, using a prime to denote a scaled item, define, for example G =

G r¯ 5 ρ 2 G

(5.52)

and in the same way the scaled matrices (G1 ) and (G2 ) . In this form (G1 ) is a function only of the aspect ratios α and β and the nondimensional spin Ω. A stress scale is given by r¯ 2 ρ 2 G, so scaled modulii are defined by K =

K r¯ 2 ρ 2 G

,

µ =

µ r¯ 2 ρ 2 G

,

κ =

κ r¯ 2 ρ 2 G

(5.53)

the entire elastic–plastic matrix by (5.47)

for the stress work terms.

C ep =

Cep r¯ 2 ρ 2 G

(5.54)

and the volume average stresses by 5.7. Total stability measures Finally, let the tensor G be the combination of the spin, gravity and stress work terms:28 G = G1 + G2

(5.48)

so that the complete test for stability, Eq. (5.36) becomes δεT Gδε > 0

(5.49)

for all admissible normal strain increments δε. If all strain increments are allowable, then the body is stable in some state if and only if the tensor G has all positive eigenvalues at that state; and if any of the eigenvalues are negative, then there is a strain increment (eignevector) that is unstable. In cases of material models with restricted strain increments, Eq. (5.49) is only tested for admissible increments. For example, if the material is incompressible, only isochoric increments are to be considered, and then the matrix can be reduced to a 2 × 2 matrix, the method is illustrated below. If it is rigid-plastic, only plastic strain increments are allowed. At a plastic state, the strain increments are either loading or unloading, with only the possibilities depicted previously in Fig. 2. 5.8. Dimensionless forms It is useful to put all of the calculations into a nondimensional form. A nondimensional spin is given as29 Ω2 =

ω2 . πρG

(5.50)

28 Note that the scalar G is the gravitational constant, the vector G the energy terms. 29 Note a factor of π that was not included in the scaling in (Holsapple, 2001), the present scaling follows Chandrasekhar.

σ¯ i . r¯ 2 ρ 2 G

σ¯ i =

(5.55)

Thus, the average stresses increase as the square of the body size. Then the scaled form of the average stresses of Eq. (5.18) are given compactly as     π σ¯ x σ¯ y σ¯ z = 5(αβ)2/3      × Ω 2 − 2Ax β 2 Ω 2 − 2Ay α 2 (−2Az ) . (5.56) Finally, the scaled stress work matrix of Eq. (5.46), involves the scaled (primed) quantities: 

G2

 ij

=

4π ¯  ep C . 3 ij

(5.57)

5.9. Relative magnitudes of terms For any state where admissible increments of strain are ¯  ep is simonly elastic, the scaled average material matrix C ij ply an incremental elastic matrix30 with terms involving scaled elastic modulii, given in Eq. (5.53). Typical values for the modulii of soils are on the order of 100 MPa. For an asteroid with average radius r¯ ≈ 10 km, mass density ρ ≈ 2 g/cm3 , the stress scaling factor is r¯ 2 ρ 2 G ≈ 30 kPa. Therefore the scaled bulk and shear modulii are on the order of 3 × 103 , while the scaled stress terms from Eq. (5.56), for any spin within the limits, are of the order of unity. Therefore, in the evaluation of the terms of the total G matrix, the ratio of the C ep elastic terms to the body force and spin terms is a large number for all asteroids with sizes up to a 30 It need not actually be isotropic, the point is only about the magnitude of the terms.

292

K.A. Holsapple / Icarus 172 (2004) 272–303

couple of hundred km diameter. Therefore, the conclusion is that, at any state allowing elastic increments of strain, and assuming typical values for the elastic modulii, the state will be stable assuming only stability as a static elastic material, Eq. (3.14) (which is certainly true for all standard models of elasticity). The stress work terms require a large amount of energy for small increments of strain, and the spin and gravity terms do not have that large energy. Only in special cases where the elastic strain energy terms are zero or negligible will stability considerations be of consequence. Of course, that is the case for fluids. Some typical applications of these concepts and results are considered next. In the detailed calculations a symbolic algebraic program such as “Mathematica” becomes essential. Since the final expressions are sometimes complex and often are found numerically, they are not all reproduced here. They come from using Eq. (5.56) for the scaled average stresses, the scaled form of Eq. (5.36) for stability, and Eqs. (5.39) and (5.57) for the evaluation of the G1 and G2 matrices, all in scaled form.

6. Maximal equilibrium configurations New and exact algebraic forms for the maximum spin for a given shape and angle of friction can now be derived. A definition of a limit state is introduced. 6.1. Equilibrium solutions There are infinitely many equilibrium configurations of a given spinning, self-gravitating body. For example, Holsapple (2001) derives general equilibrium stress distributions for such ellipsoidal bodies, and obtains a particular solution and three linearly independent homogeneous solutions. That is a consequence of the fact that there are three equilibrium equations for six stress components. Any linear combination of the three homogeneous solutions added to any particular equilibrium solution gives another equilibrium solution. 6.2. Elastic equilibrium solutions Linear elastic analyses of this problem (e.g., Dobrovolskis, 1982, and others) calculate a one specific equilibrium state by (implicitly) assuming that the body with no gravity or spin forces is in a state with zero (residual) stresses. That assumption fixes the contribution of the homogeneous solutions. Then they compare those elastic stresses to some yield criteria, to find a condition for first yield somewhere in the body.31 Those elastic stress states are not valid if there 31 In fact, for the Mohr–Coulomb yield criteria with no cohesion, the linear elastic analyses gives that the criteria is violated at the equator for even zero spin, indicating that the body must deform plastically as it accretes, even if not spinning. Or, one is forced to assume cohesion at the surface.

are any residual stresses. That limits the usefulness of elastic analyses, whether linear or not. 6.3. Plastic limit equilibrium states Here the limit configurations are found where all points of the body are at the plastic yield criteria. No larger spin can exist without global body disruption or shape readjustments. That is a limit-state concept as used in limit analysis of plasticity theories.32 It will be shown that those states are unique, are independent of any elastic behavior and existing residual stresses, and determine precise spin limits of the body for a given shape and density. Suppose then that all points in the body are just at the plastic limit.33 Then at all points, the maximum and minimum principal stresses satisfy the yield criteria Eq. (4.5); so, integrating over the body, the average maximum and minimum principle stresses satisfy the same criteria. The complete stress distributions as a function of x, y, and z for this problem were given in (Holsapple, 2001, Eq. (10)). An important feature was that the shear stresses were identically zero, and the three normal stresses σx , σy , σz are in the same ratios throughout the body. As a consequence, those three are also principle stresses, and one is maximal throughout the body, and another is minimal. Therefore, for a given shape and spin, the average of the maximum principal stress is also the maximum of the average of a particular normal stresses.34 However, for different spins and shape states, the maximum can be any one of the three and the minimum one of the other two. That gives six possibilities to consider. All six possibilities can be considered together if a set of three constants λx λy λz  is defined by the rule  m if x is the direction of maximum principal    stress,   0 if x is the direction of intermediate λx = (6.1) principal stress,       −1 if x is the direction of minimum principal stress, with similar rules for the y and z components. For example, suppose that the average σ¯ x is the maximum (least compressive) and that σ¯ y is the minimum (most 32 What are found here are equilibrium states that are everywhere at or below yield. Then the spin loads for collapse are no less than those present at such a state. To prove they are also the maximum, it is necessary to demonstrate a kinematically possible velocity field, which is presented in (Holsapple, 2001). 33 The fact that for the cohesionless Mohr–Coulomb body there is a limit state with the stresses at all points just at the yield was proved in (Holsapple, 2001), where the complete stress distributions were derived. This is not a property to be expected for other material models, and is not a general feature for limit analysis. In other cases the maximum possible spin could occur with just some points in the body at yield. 34 If the maximal principal stress were a different stress component at different points in the body, then the average of the maximal principal stress would not be the average of a particular component.

Equilibrium figures of spinning solid bodies

compressive). In that case λx λy λz  = m −1 0. In cases where two or more stresses are equal, then Eq. (6.1) gives more than one possibility to consider. Using those parameters, all global yield states defined by averaging Eq. (4.5) satisfy the same simple form λx σ¯ x + λy σ¯ y + λz σ¯ z = 0,

(6.2)

where one of the λi terms is zero, one is −1, and one is equal to m. Using the average stresses from Eq. (5.56), this can be explicitly solved to give the spin at that global yield state; the simple algebraic result is Ω2 =

2(λx Ax + λy y + λz λx + β 2 λy β 2A

α2 A

z)

.

(6.3)

This Eq. (6.3) is the general criteria that define the limit yield states, and is a major result of this paper. It gives the maximum possible spin for which all points in the body are at yield. The right-hand side of this equation is some definite function of the two aspect ratios α and β, and of the single material parameter m. No considerations of residual stresses or of elastic behavior within yield were needed. The material may be compressible or incompressible. No plastic flow rule is needed. The only complexity is that the three λx , λy , and λz , according to Eq. (6.1) differ according to the relative magnitudes of the principle stresses, which are different in different regions in the α, β plane. The curves that delineate those regions can easily be found by equating two of the average stresses from Eq. (5.56) and solving for the common scaled spin. For the convenience of the reader, the six possible cases are listed in Table 1. Note that all of the results for limit spins are size-scale invariant: i.e., they are the same for all body sizes. Results from these equations are now presented.

293

7.1. Limit spins for fluid bodies For a fluid, the angle of friction is zero, so that the parameter m = 1. At each point in the body the three normal stresses are equal (the negative of the pressure). Then Eq. (6.3) with all six cases can be used, but since Eq. (6.3) is invariant to the overall sign of the λi , only three are distinct. Therefore apply the cases 4, 5, and 6 of the Table 1:  1 −1 0, λx λy λz  =

0 1 −1, 1 0 −1,

(7.1)

to get three equations for the spin: Ω2 =

2(Ax − β 2 Ay ) , 1 − β2

Ω2 =

2(β 2 Ay − α 2 Az ) , β2

2(Ax − α 2 Az ) . 1 The last of these gives that   Ω 2 = 2 Ax − α 2 Az Ω2 =

It is instructive to begin with the fluid cases, which are included in this present analysis by taking the angle of friction equal to zero, and to compare the results to the classical results. Both the possible spins and their stability are considered. Table 1 Case

Stresses in the body satisfy

The λi

1

σx < σy < σz

−1 0 m

2

σx < σz < σy

−1 m 0

3

σy < σx < σz

0 −1 m

4

σy < σz < σx

m −1 0

5

σz < σx < σy

0 m −1

6

σz < σy < σx

m 0 −1

Formula for Ω 2   2 Ax − mα 2 Az  Ax −mβ 2 Ay  2 mβ 2 −1   2 2 Ay − mα2 Az β  mAx −β 2 Ay  2 m−β 2   2 2 Ay − α 2 Az mβ   2 2 Ax − αm Az

(7.3)

and the other two are not independent. Either can be used to give a consistency constraint   β 2 (Ax − Ay ) + α 2 1 − β 2 Az = 0. (7.4) This Eq. (7.4) is a required relation between the two aspect ratios α and β, otherwise equilibrium is not possible. The ellipsoids satisfying that relation are called the “Jacobi ellipsoids.” For any given α, one could (numerically) solve for the corresponding β that satisfies that constraint. The combinations of aspect ratios of the Jacobi ellipsoids are shown here in Fig. 5. Chandrasekhar (1969) has the consistency relation Eq. (7.4) in terms of a quantity Axy as a 2 b2 Axy = c2 Az

7. Fluid bodies

(7.2)

(7.5)

a form he attributes to Jacobi. It is easily shown that the Axy is given as Axy =

1 a 2 (1 − β 2 )

(Ay − Ax )

(7.6)

so Eq. (7.5) is the same as Eq. (7.4). Chandrasekhar gives on his p. 103 a table for the consistent α and β pairs that satisfy Eq. (7.4), the same as shown on Fig. 5 here. Then, for any particular shape combination, one can determine (numerically) the equilibrium spin from Eq. (7.3); those are also depicted in Fig. 5. Equations (7.3) and (7.4) are, according to Chandrasekhar (1969), almost exactly the same forms given originally by Jacobi. Various integral forms for the spin are possible using both of Eqs. (7.3) and (7.4) and Eq. (5.4). One such form is ∞ Ω = 2αβ 2

0

u du (u + 1)3/2 (u + β 2 )3/2 (u + α 2 )1/2

(7.7)

294

K.A. Holsapple / Icarus 172 (2004) 272–303

Fig. 5. The shapes and spins of the Jacobi ellipsoid fluid bodies. For each β = b/a, only a single α = c/a is permissible, and only a single spin is in equilibrium. These results date to Jacobi in 1834.

which is a nondimensional version of the form given for the Jacobi ellipsoids by Chandrasekhar (1969) (p. 101, Eq. (5)), using his definition of B12 as this integral. While the integral Eq. (7.7) requires elliptic functions for evaluation, a program such as “Mathematica” can easily numerically determine the β value for any given α from Eq. (7.4) and then evaluate (7.7) to get the spin for that α. A special case of the results just given are for an oblate spheroid with a = b > c, and more special results are possible. Then Ax = Ay and β = 1, and the geometry constraint Eq. (7.4) is identically satisfied, so that there is a solution for every α. Using Eq. (5.3) in Eq. (7.3) gives the spin as   Ω 2 = 2 Ax − α 2 Az = 2Ax − 2α 2 (2 − 2Ax )   = 2 1 + 2α 2 Ax − 4α 2 . (7.8) Further, in this case the gravitational integral has a relatively simple form: αCos−1 (α) α2 − (1 − α 2 )3/2 1 − α 2 and a little algebra gives the equilibrium spin rate as

Ax =

Ω2 =

2α(1 + 2α 2 ) 6α 2 Cos−1 (α) − . 2 3/2 (1 − α ) 1 − α2

(7.9)

(7.10)

√ When rewritten in terms of the ellipticity e = 1 − α 2 this is exactly the formula given by Maclaurin for spin rates of these oblate fluid bodies, known as the Maclaurin spheroids. These results will be plotted below with the results for a solid body. 7.2. Stability of fluid shapes The limit curves for fluids are unique equilibrium states with one spin for each shape. An important question is of

the stability of these states. Suppose there is some event that perturbs the state. For example, suppose there is an impact event. Of course, since we only consider ellipsoidal shapes, we cannot consider the details of cratering or other configuration modifications, but we can consider the stability to a small change in the overall ellipsoidal shape. Such perturbations give necessary conditions for stability, not sufficient. The presentation of stability was given above in Section 4 using a perturbation with constant strains. For incompressible fluids the strain energy terms are identically zero,35 and stability involves only the G1 (or scaled G1 ) matrix defined in Eq. (5.39), and for stability, it is necessary that δεT G1 δε > 0

(7.11)

but restricted to permissible strain increments. Since the fluid materials are assumed to be incompressible, the strain increments must be isochoric. Then the three normal strains can be expressed in terms of two arbitrary x and y strains as    1 0 γx = Qδ εˆ δε =  0 1  (7.12) γy −1 −1 with the matrix   1 0 Q= 0 1 . −1 −1

(7.13)

The δ εˆ strain matrix has only the two independent components γx and γy . Using Eq. (7.12) in Eq. (7.11), the state will be stable in a state if the 2 × 2 reduced matrix G3 = QT G1 Q has positive eigenvalues at that state. 35 The stress work increment of an inviscid fluid is just −pδv, and for incompressibility the volume change is zero.

Equilibrium figures of spinning solid bodies

295

Fig. 6. The eigenvalues of the Maclaurin spheroids and the Jacobi ellipsoids. The Maclaurin spheroids have β = 1.0. For the Maclaurin spheroids one eigenvalue is positive everywhere, but the other is negative below α = 0.5827. There the Maclaurin spheroids are unstable with the unstable (strain) eigenmode 1 −1 0. For each β the Jacobi branch has only the aspect ratio α shown in Fig. 5. Those are all stable for the constant strain increments considered here. The configuration at β = 1 and α = 0.5827 has a zero eigenvalue and is neutrally stable; it is the point of bifurcation of the Jacobi ellipsoids from the Maclaurin spheroids.

Therefore, to test stability at a given limit state, with given aspect ratios α and β, the eigenvalues of the reduced matrix G3 (evaluated with the limit spin) are investigated, and a negative value indicates instability. As the first example case, consider the Maclaurin spheroids, for which β = 1. Using Eq. (7.10) for the spin, the two eigenvalues of the G3 matrix have been found (using “Mathematica”) and are plotted as a function of the aspect ratio α in Fig. 6. These results are very interesting in comparison to the stability results reported by Chandrasekhar (1969). From his spectral stability analysis he finds that the fluid bodies are unstable below α = 0.3303, neutrally stable from that value to α = 0.5827, and stable for α greater than that. The neutral and unstable modes had (oscillating) strains of the same eigenmode 1 −1 0 found here. For his neutral stability range, he shows that any fluid viscosity makes those modes unstable, so they are states of secular instability. Of course it is known that spectral stability does not imply other definitions of stability. With the definition of stability used here, all of the configurations with α < 0.5827 are unstable, irrespective of considerations of viscosity. The additional points of bifurcation Chandrasekhar finds for the third harmonics are of no interest here since they occur for aspect ratios below the present stability limit. The Jacobi ellipsoids were also analyzed here, the results are also on Fig. 6. For them the smallest eigenvalue has a zero value at the state β = 1.0, α = 0.5827, the point at which the Jacobi ellipsoids bifurcate from the Maclaurin

curve, and are positive below that. See Chandrasekhar for a complete discussion of this bifurcation. Then, the Jacobi ellipsoids are stable by the present analysis. For the Jacobi ellipsoids, Chandrasekhar also finds the neutral point where they bifurcate from the Maclaurin spheroids at α = 0.5827. For second-order oscillations, (linear displacements, the same as considered here) he shows that they are all stable below that point. However, he goes on to consider third-order variations (quadratic displacements), and re-derives the result of Poincaré. There is a bifurcation to a Pear-shaped body at the shape α = 0.345 and β = 0.4322, but Chandrasekhar finds that these shapes are unstable. Considerations of higher-order perturbations and possible bifurcations to pear-shaped (and perhaps binary?) bodies by the present approach are left for future study. However, it is noted that for compressible solid bodies, higher-order displacement increments will lead to nonuniform strains and bodies that are not at uniform density, so the calculation of the gravity potential increment becomes problematic.

8. Limit states and stability for solid bodies I now return to the solid case,36 and use Eq. (6.3) with a nonzero friction angle. The additional complexity then is that various combinations of α and β will have different 36 Again recall that it is assumed that the cohesion is zero, an assumption appropriate for rubble-pile bodies. The additional complexities of including cohesion will be addressed in a subsequent paper.

296

K.A. Holsapple / Icarus 172 (2004) 272–303

combinations of the x, y, or z stress components as maximum and minimum, and generally each case must be considered in turn.37 The spin is always around the shorter, z-axis. Clearly for small spin rates gravity forces dominate and the x diameter, being longer, has more compression than the other two. Thus, the smallest spin conditions always have σx the smallest (most negative) and σz the largest stress, which is case 1 of Table 1. Then for larger spins, the stress along the xand y-axes become less compressive due to the centrifugal forces, and other cases of Table 1 are possible. Finally, for the largest possible spin rates the smallest stress (most compressive) will ultimately be σz along the z-axis, giving one of cases 5 or 6.

For oblate bodies a = b, β = 1.0, and σx and σy are equal. First consider small spin rates; then these two stresses are the smallest and σz is the largest principle stress. In this case, Eq. (6.1) has two possibilities  0 −1 m, λx λy λz  = (8.1)  −1 0 m.

The curves from Eq. (8.3) can be shown on a plot of scaled spin versus aspect ratio. In (Holsapple, 2001) these curves were presented, although the analysis was quite different and less direct. There the entire stress state was first determined. Then the yield criterion was used to get a required friction angle as a function of the spin and aspect ratios. The curves for a constant friction angle were then obtained by contour plots of constant friction angle values. Here the general explicit algebraic formula of Eq. (6.3) has been discovered, and in this case of the lower region (small spins) for an oblate spheroid, it gives the explicit results of Eq. (8.3). Next larger spin rates are considered. At some spin, the equal x and y normal stresses also equal the z stress, which is exactly the fluid case discussed above. Therefore the result for the Maclaurin spheroids for fluid bodies is also the curve demarcating the upper and lower regions of possible spin of a solid body. Then above the Maclaurin spheroids cases (larger spins), the x and y components become the largest, and the z the smallest (most compressive). Then, using either of  m 0 −1 , λx λy λz  = (8.6) 0 m −1 ,

Since β = 1, and Ax = Ay , both choices give the same result, the case 1 of Table 1, that   Ω 2 = 2 Ax − mα 2 Az . (8.2)

the cases 5 and 6 of Table 1, gives the simple formula:   α2 Ω 2 = 2 Ax − Az . (8.7) m

This is the general result for a solid when σx is the smallest stress and σz is the largest. Using Eqs. (5.3) and (7.9) gives the explicit solution for the limit spin rate as

This is the general result, whether or not the body is oblate, for any case when σx is the largest stress and σz is the smallest. Using Eqs. (5.3) and (7.9) now gives

8.1. Limit spin of solid oblate spheroids

Ω2 =

2α(1 + 2mα 2 ) 2(1 + 2m)α 2 −1 Cos (α) − . (1 − α 2 )3/2 1 − α2

(8.3)

This reduces to the Maclaurin formula when m = 1. Therefore this is a generalization of the Maclaurin formula to cohesionless solids. It is a lower branch of spin limits only and is below the Maclaurin curve on a plot of spin versus aspect ratio. It can easily be evaluated for any value of the material parameter m and the aspect ratio α. Further, it is obvious that, unless 2α(1 + 2mα 2 ) 2(1 + 2m)α 2 Cos−1 (α) > 2 3/2 (1 − α ) 1 − α2 there is no solution, and when √ (1 + 2m)α 1 − α 2 −1 Cos α = (1 + 2mα 2 )

(8.4)

(8.5)

the spin is zero. Thus there is a limit to the oblateness for a solid body without spin, depending on the friction angle. 37 I again note that the use of the Drucker–Prager yield criteria gives nearly identical results, and does not have this complexity. Instead the algebra is more complex, but the use of “Mathematica” overcomes that difficulty.

Ω2 =

2α(m + 2α 2 ) 2(m + 2)α 2 −1 . Cos (α) − m(1 − α 2 )3/2 m(1 − α 2 )

(8.8)

This is a generalization of the Maclaurin result for the limits of the larger spins for solid bodies above the Maclaurin curve. In this case, there are possible spins for all α values between 0 and 1. The plots of these results are shown in Fig. 7, which also shows the previous Maclaurin curve for a fluid from Eq (7.10). The plots of Eq. (8.3) are the curves in the lower regime of that figure, beneath the curve for the Maclaurin spheroids. The curves above the Maclaurin curve are given by Eq. (8.8). The stability curve shown on this figure will be discussed below. 8.2. Other ellipsoidal shapes The general spin limit Eq. (6.3) can be applied to any ellipsoidal shape to find the maximal spin rates. The potential terms Ai may require numeric evaluation. However, for a prolate body closed-form expressions are again available, but now using log functions. However those are not given here. For prolate solid bodies, there are three orderings of the stresses to consider. The results are shown in Fig. 8. These

Equilibrium figures of spinning solid bodies

297

Fig. 7. Spin curves and stability limits for an oblate solid with various angles of friction. States on the limit curves are at the maximum spin possible for that angle of friction, and the entire body is just at the Mohr–Coulomb yield limit. For each friction angle, there are two curves, one upper limit and one lower limit and all states between are possible equilibrium states. The Maclaurin curve in the middle is along shape combinations where all stresses are equal, as for a fluid. The stability limit curve shown assumes a loss of shear rigidity for small strain increments and reduces the allowable shapes considerably.

Fig. 8. The equilibrium limits for a prolate solid body with various angles of friction. For each angle there are two curves, an upper limit and a lower limit. Those are from the Eq. (6.3), recognizing the different ordering of the stresses in the three regions. The stability curve is discussed below.

are √ just as in (Holsapple, 2001), but scaled differently by a π. What are new here are the specific formulas for those curves. For any aspect ratio values α and β, the spin limits can always be determined in a definite way. First curves demarcating the regions of stress ordering can be found by equating the relations for the scaled average stresses in Eq. (5.56). Then in each region, the appropriate case from Table 1 can be found. As a final √ example, I assume that the intermediate axis b = ac. This case has all six stress combinations present. The resulting equilibrium limit curves are shown in Fig. 9.

8.3. Stability of solid bodies A significant property of the incompressible inviscid fluids is that the strain energy term is identically zero, which greatly simplifies the stability analysis. However, for a solid body, the stress work terms must be included. At the plastic limit state when loading, those are from plastic strain increments. Or, for states inside the plastic limit curves or for unloading, the increments of strain are elastic. With representative modulii values it has already been noted that the magnitude of those elastic terms is positive and generally they dominate the magnitudes of the spin and gravity terms. As a consequence, any state with elastic changes is stable

298

K.A. Holsapple / Icarus 172 (2004) 272–303

√ Fig. 9. Equilibrium limit curves when the body has the “potato” shape described by b = ac. The three darker curves demarcate the six regions of different stress value combinations. Their intersection at the point with α = 0.473 is at a Jacobi ellipsoid for a fluid where all three stresses are equal. The thick shaded curves are stability limits when the shear stiffness goes to zero, as discussed below.

for representative values of the elastic modulii. For a solid body that is generally the case for all states within the limit spin curves. Exceptions occur only when the stress work terms are small, zero or negative. Those terms are a product of stress increments times strain increments integrated over the body. Constitutive behavior gives the relation between the strain and stress increments. They are zero in two limits: either the stress increment is zero for some strain increments, or the strain increment is zero for some stress increments. The first case requires modulii (e.g., shear and bulk modulii) to be zero; an example is for the shear of an inviscid fluid where any shear strain has zero shear stress. In a solid one might imagine a loss of elastic stiffness temporarily during an impact event or a tidal encounter, or perhaps by some long-term creep mechanism. Since nothing is known of the actual properties of asteroids, it is appropriate to consider all possibilities. In the second case, the modulii become infinite, so that the material becomes rigid, for at least some stress increments. A solid material example is rigid-plastic behavior, for which elastic strains are zero for increments of stress within yield, or when at yield but unloading. However, first some important points must be made. Here only constant strain increments over the body have been introduced. It is well known in plasticity theories that bodies that are everywhere at yield often do not have additional plastic strains that are uniform; but instead, small variations in properties or state cause plastic strains to “localize” along surfaces, leading to “shear banding.” That can happen when the work hardening is small or negative. These phenomena cannot be included in the present analysis, and certainly need further study. The grooves seen in many Solar System bodies are highly suggestive of such processes.

Also, any state within the limit spin curves can be expected to have both regions at the plastic yield limit, and regions within yield. Therefore, the assumption of a rigidplastic material that the elastic strains are zero in selected regions cannot be accommodated within the constant strain assumption. That would require nonhomogeneous strains. The assumption that a shear stiffness is zero only in selected regions cannot be easily studied. One cannot study loss of shear stiffness in selected regions, for example along a shock surface in an impact. One can only make the strong assumption that it is zero everywhere. Therefore, the following study of certain limiting cases is only scratching the surface of possible stability considerations for a solid body. However, it may be indicative of the extent of states where stability considerations are important. So I consider exceptions, when the stress work terms are small, zero or negative. That includes plastic increments on the limit curves; and when within those limit curves, includes elastic increments with the elastic modulii zero or near zero, so there is little or no elastic stiffness. Recall that it was mentioned above that all analyses here use the incremental stress-strain relations Eq. (4.12) specialized for a general Mohr–Coulomb plasticity model with work hardening and with incremental elastic unloading. That gives many possibilities of combinations of the elastic constants K and µ (and others for an anisotropic material), the angle of friction parameter m, the flow rule determined by the parameter m1 , and the work-hardening (softening if negative) modulus function κ. All of these affect the incremental elastic–plastic stress-strain relation, and thus the stress work increment. Here I shall investigate a few combinations that seem to be of the most interest and lead to important stability considerations.

Equilibrium figures of spinning solid bodies

299

Fig. 10. The scaled stability measure for plastic strain increments and states on the plastic limit curves with a friction angle of 30◦ .

8.3.1. Stability at limit spin states, oblate solid bodies I shall consider in turn several possibilities of material parameters. First is for an oblate body. Suppose a limit spin state, so the body is entirely at the plastic limit, and assume a plastic increment of strain according to the flow rule, the outward vector shown in Fig. 2. Then the state must remain on the yield function, and, if the body is perfectly plastic, the stress increment is zero, so the elastic strain increment is also zero. This suggests an interesting case for a stability test. However, the plastic strain increment, which must be outward to the yield function, is of the form38 m1 0 −1. While the stress work increment is then indeed zero for those strain increments, the unstable eigenvector of the G1 matrix by itself (involving the spin and gravity terms) is 1 −1 0 and not in the same direction as the plastic strain increment. So the fact that the stress work increment is zero does not imply a stability limit. Instead, in this and other cases, the entire G matrix with combined spin, gravity and stress work terms must be investigated for a common assumed strain increment. So, I consider a state on a limit curve, a plastic strain increment, and include the possibility of work-hardening (or softening). A plastic loading increment must be something like δε = λm1 0 −1

(8.9)

for some positive λ, although the location of the entries depends on the relative stress magnitudes. Then the value of the stability measure δεT G δε (the scaled version of Eq. (112)), evaluated with a limit spin given by Eq. (6.3) in the appropri38 In this oblate case with a = b, the x and y stresses are equal, so the state is actually at a corner of the Mohr–Coulomb yield function. That admits two possible directions, the other is 0 m − 1, and also some vector between these two. However, in all cases the important results are found with just one of these directions; it is comforting that the singular corner behavior of the yield function is of no consequence.

ate stress regime is some function of the aspect ratio α times λ2 . For stability the function must be positive. As an example, for oblate bodies in the region above the Maclaurin curve the limit spin is given by Eq. (8.7). I take the friction angle to be 30◦ , with m = 3 and m1 = 3, the scaled bulk modulus to be 104 , and first, the perfectly plastic case of κ = 0, zero hardening. The stress work terms are zero with this strain increment irrespective of the elastic constants, since the elastic strain increment is zero and, as previously noted in Section 3, the plastic stress work is zero in this perfectly plastic case. Then the value of the stability measure is that of the spin and gravity only and is as shown in Fig. 10. It is positive for all α, except those below 0.150937, so all except the most elongated shapes are stable for a perfectly plastic increment of strain on the limit curves when the friction angle is 30◦ . However, if one allows for strain softening, a very different result emerges. The scaled stability measure with only the spin and gravity terms have the values of 20 or less in this figure. But when the scaled hardening (or softening) parameter κ  is included, the (scaled) stress work increment terms are given in this numerical example by the relation  −1 4π  3κ  κ 1+ . (8.10) 3 (12K  + 52µ ) If the work-hardening parameter is negative, the stress work terms release energy upon strain increments. For the stress work term to release enough energy to offset stability, it only need be negative of the order of −10 to −20, which using Eq. (8.10) happens when the scaled hardening modulus is negative of the order of −2 to −5. On Fig. 10 is plotted the result of the value of the stability measure using scaled softening of κ¯ = −2 and of κ¯ = −4.7. The first value makes the solid oblate body unstable below about α = 0.35, while the latter is sufficient to make all oblate bodies unstable at the plastic limit. The unstable mode is the plastic flow vector δe = m 0 −1. That mode

300

K.A. Holsapple / Icarus 172 (2004) 272–303

is a lengthening of the x-axis with, assuming φ = 30◦ , 1/3 as much shortening of the z-axis, with an associated volume bulking of 2/3 of the percent of x lengthening. The use of a nonassociated flow rule, with m1 = 1 gives no volume change but does not affect the general nature of this result, but only the numerical values by a small amount. With a stress scale factor of about 30 kPa for a 10 km diameter asteroid, the required softening modulus value needed to completely upset stability is about −150 kPa, compared to the value of shear modulus mentioned above of 40 MPa. Therefore, to completely offset stability at the limit spins requires a softening modulus of only 0.4% of the shear modulus, much less than the couple of percent noted for the experimental data discussed above. While the elastic modulii are large compared to the scaled spin and gravity terms, to affect stability the work-hardening modulus need only be negative on the order of the spin and gravity terms, or a very small part of the elastic constants. Therefore, the conclusion is that all states on the limit curves become unstable with only a very small amount of strain softening, which is a standard feature of soil-like materials. The unstable mode would be an elongation of the largest axis, with a shortening of the spin axis.39 That assumes a uniform strain. As mentioned above, the more probable result would be a strain localization and shear banding along planes inclined to the longest axis by a steep angle (associated flow rule) or by as little as 45◦ (nonassociated flow rule, no volume change). 8.3.2. Loss of elastic stiffness, oblate bodies What about stability for states below the limit spin states? The results above predict instability on the limit curves when there is strain softening, but says nothing about states inside the limit curves. One rather extreme assumption discussed above is that the elastic modulii are zero. So, I now suppose there is some unspecified mechanism makes the elastic modulii become zero. Whether this can actually happen is unknown, but it is interesting to contemplate and study.40 Specifically, I discard the stress work increment terms, and study the consequences regarding stability. The approach is just like that for the fluid, but instead of using the limit spin relations for a fluid, those for a plastic solid are used. Basically, the eigenvalues of the matrix G1 are investigated for sign, now using the spin relations for a solid. For 39 Does the body disrupt, or does it go to some new but stable configuration by a large global change of configuration? The actual deformation event cannot be predicted from the present infinitesimal stability theory, but one might wonder about the formation of binary bodies. 40 In cratering mechanics, some like to invoke a mechanism called “acoustic fluidization” which is a late-stage loss of stiffness replaced by a viscous mechanism. That ad-hoc assumption can produce late-stage crater modification and slumping. However, that mechanism requires loss of strength long after the passage of the shock due to the impact, which is much harder to imagine than loss during the shock passage. But here, vibrations in the shock due to a noncatastrophic impact might be a mechanism for a temporary loss of stiffness.

oblate bodies in the region above the Maclaurin curve, that spin relation is given by Eq. (8.7), and for states below the Maclaurin curve is given in Eq. (8.2). The G1 matrix has three eigenvalues, which are complicated functions of the body aspect ratio α and the scaled spin Ω. The eigenvectors are a 3 × 1 matrix of normal strains. It is found by numerical evaluations that one of the eigenvalues is negative in all cases, with an eigenvector being a strain vector with a volume increment. Therefore, if there is no stiffness in bulk, the body is simply unstable in all states and will collapse. However, as long as a bulk stiffness exists that is of the same order of the stresses or higher, that mode becomes stable. The limiting case of incompressibility, having infinite bulk modulus, does not allow that unstable mode. That case, with either finite or infinite bulk stiffness is then studied further. When the strains are restricted to incompressible modes, using the transformation matrix of Eq. (7.13) the G1 matrix can be reduced to the G3 matrix with only two eigenvalues. The corresponding physical assumption is sufficient elastic stiffness in bulk, and zero stiffness in elastic shear. For the case of a fluid with the yield parameter m = 1 and assuming incompressibility, the stability limit on the Maclaurin curve was previously found at α = 0.5827. For cases of solids, with m > 1, one can solve for the maximum spin for which both the two remaining eigenvalues are positive as a function of the aspect ratio α. The limit spin for each aspect ratio α have been determined in that way (using “Mathematica”), and overlaid on the previous plots of limit spins in Fig. 7 above. These curves do cut across the region of interest. For each aspect ratio there is then a maximum spin for stability, and for greater spin the body is unstable. The mode of the instability is again the strain 1 −1 0. This is a shear strain at 45◦ of the principal x- and y-axes. Therefore, only the shear modulus must be zero for this instability; the bulk modulus plays no role. All points above the stability limit curve are unstable if the shear modulus is zero. Limiting curves for nonzero but small shear stiffness have been found also, but whenever the shear modulus is on the order of the stress or larger, all states within the plastic limit state become stable. Thus those intermediate cases are not considered further. It is seen that loss of shear stiffness can substantially reduce the possible spins. For a spherical oblate body the maximum scaled spin becomes Ω = 0.73, which for an asteroid with mass density of 2 g/cm3 is a period of 3.7 hr/cycle, considerably slower than the “spin limit” of 2.1 hr/cycle often quoted. Elongated bodies have this stability spin limit reduced to even slower values as shown on Fig. 7. 8.3.3. Loss of elastic shear stiffness, prolate bodies Since most asteroids have an elongated shape, the prolate case with two short axes and one long one is usually of the greatest interest. In this case, I will ignore the first two cases discussed above for an oblate body, and move directly to the case where there is a loss of shear stiffness. Again the con-

Equilibrium figures of spinning solid bodies

301

Fig. 11. Scaled spins of over 1000 asteroids with diameter > 2 km versus their aspect ratio α = c/a, assuming prolate ellipsoidal shapes. Most are beneath the stability curve assuming no shear stiffness, but there are a noticeable number above. Those asteroids would have to have shear stiffness or cohesion to remain in those states. They are all relatively small S-types or “other” types, but none are C- or M-types.

siderations are for the values of the eigenvalues of the G1 matrix. Since again loss of bulk stiffness always makes the body unstable, the case of zero shear stiffness is the one of interest. Again “Mathematica” is exercised to find the values of the remaining two eigenvalues of the G3 matrix as a function of spin and aspect ratio. For the more spherical bodies, with an aspect ratio α of about 0.5 or greater, there is an upper limit on spin for stability. That upper limit curve is shown on the previous Fig. 8. For the spherical case, that limit is again Ω = 0.73 while for slightly elongated bodies a value of about 0.8 is permissible. For very elongated bodies, there is also a lower limit on spin as shown, but that case is of little practical interest. 8.3.4. Comparison to asteroid spins The results above can be compared now to the actual data for asteroid spins and shapes. Generally, the observed light curves allow the estimate of the ratio β = b/a of the intermediate and longest axis, but not easily the aspect ratio α. Therefore, one must make some assumption about the relation between α and β in order to plot the data. I assume now that all asteroids are prolate, so that α = β. I also must make some assumptions about mass density in order to scale the raw data. I assume values, respectively for the C-types, S-types, M-types, and all others of 1.5, 3.0, 8.0 and 3.0. Then the data41 are plotted over the curves determined by 41 The data is from the database at URL http://www.asu.cas.cz/ ~asteroid/lcsumdat20031219.txt, compiled by A. Haris (P. Pravec, private communication).

the theory, as shown in Fig. 11. These data do not include the recently discovered “fast-rotators,” which would plot off the top of this plot, and clearly have to have some cohesion. A subsequent study will address that question. If the smallest axis a of any one were less than the intermediate one b, then the data point would displace to the left. If the mass density were smaller than assumed, then the data point would be displaced upward. Most asteroids are beneath the stability curve assuming no shear stiffness, but there are about 40 above that curve. Those asteroids would have to have shear stiffness or cohesion to remain in those states. They are either S-types or “other” types, but none are C-types or M-types.

9. Concluding comments Significant progress has been made on the problem of the equilibrium and stability of spinning, self-gravitating solid bodies. The maximal spins where a cohesionless Mohr– Coulomb solid ellipsoidal body is in equilibrium with all points at the yield limit have been found from a simple closed-form algebraic formula, Eq. (6.3). That formula was found using a general result for the average stresses in a spinning, self-gravitation body of any material, Eqs. (5.18) and (5.56). The case of fluid bodies is included in this equilibrium analysis. The stability question of these solid bodies is complex. Constitutive equations for elastic–plastic materials are not smooth at yield, so classical stability analyses such as the

302

K.A. Holsapple / Icarus 172 (2004) 272–303



 spectral method, which relies on linearization for small oscillating increments from an equilibrium state, are ruled out. Considerable effort was devoted to developing a new measure of stability. A new approach using a measure depending only on the equilibrium state and defined displacement increments has been presented. Its application to other definitions of stability for special materials is given to show that in those cases, it is consistent. However, there are still unknown features about that measure. In the application of this approach to spinning Maclaurin spheroid fluid bodies, I recover the region above α = 0.5827 (eccentricity e = 0.8127) of stability, but obtain instability below that value. That is in contrast to classical results which have stability down to α = 0.3033 (e = 0.9529), unless there is a viscous shear stress. Whether this difference is due to the fact that the present approach has fluids as a limit of a solid that admits shear stresses, or is a difference in the stability criteria is not yet known, although I favor the first reason. For general solid bodies, I discover two cases of instability. On the limit curves, the bodies are plastically unstable for very small strain softening, and the predicted mode is a lengthening of the body along its longest axis. In fact, the theory only admits uniform straining, but the failure would most probably be by strain localization along inclined planes. States inside the limit curves for a given friction angle are stable with conventional values of shear modulus. However, if the bulk modulus is zero, all such states are unstable; and if the shear modulus became zero, there is a definite curve of spin above which the states are unstable. Those curves are depicted in Figs. (7), (8), and (9) for different shape ellipsoids, and the data on asteroids is compared to those curves in Fig. (11). Future research is underway to investigate the importance of cohesion, the equilibrium and stability when there are tidal stresses, and further clarify stability questions.

Appendix A. Forms of the stability measure Equation (3.7) of the paper gives the condition for stability using integrals of variations    R T δT : δH dV − δb · δu dm + δas · δu dm  > δtR · δu dA. (A.1) Other useful forms of this expression are possible. Consider the following:    δ TR : δHT dV − δ b · δu dm + δ a · δu dm  − δ tR · δu dA    R T R 2 T = δT : δH dV + T : δ H dV − δb · δu dm

− −



b · δ 2 u dm + δtR · δu dA −



 δa · δu dm +

a · δ 2 u dm

tR · δ 2 u dA.

(A.2)

Since the theorem of virtual work can be applied to any field, all of the δ 2 terms cancel so    R T δ T : δH dV − δ b · δu dm + δ a · δu dm  − δ tR · δu dA    = δTR : δHT dV − δb · δu dm + δa · δu dm  − δtR · δu dA. (A.3) Now decompose the acceleration terms and use    δ aˆ · δu dm = δ aˆ · δu dm + aˆ · δ 2 u dm  = δ aˆ · δu dm

(A.4)

since aˆ is zero in the initial configuration. So the aˆ terms in Eq. (A.3) cancel and Eq. (A.1) gives that    δ TR : δHT dV − δ b · δu dm + δ as · δu dm  > δ tR · δu dA. (A.5) Therefore, the definition of stability using the Piola stresses can have the variation symbols either inside or outside the integrals. Next, forms with the Cauchy stress are derived. The stress work terms can be directly written with the Cauchy stress:   R T T : δH dV = T : δe dv,   tR · δu dA = t · δu da, (A.6) so that Eq. (A.5) also gives    δ T : δe dv − δ b · δu dm + δ as · δu dm   δ t · δu da.

(A.7)

When the surface terms are zero, either of Eqs. (A.5) or (A.7) are δ 2 SW + δ 2 PE + δ 2 KEss  0.

(A.8)

Moving the variation inside the stress work term of Eq. (A.7) gives three terms:  δ 2 SW = δ T : δe dv    = δT : δe dv + T : δ(δe) dv + Tr[δe]T : δe dv (A.9)

Equilibrium figures of spinning solid bodies

using δ(dv) = δJ dv ≈ Tr[δe] dv. The term δ(δe) = δ × (grad δu) requires special attention, since the δ operator does not commute with a spatial gradient. Using the present as reference, then   δ(δe) = δ(grad δu) = δ (Grad δu)F−1   = Grad δ 2 u F−1 + (Grad δu)δF−1     = grad δ 2 u − (Grad δu) F−1 δFF−1   = grad δ 2 u − grad δu δFF−1 = grad δ 2 u − grad δu grad δu = grad δ 2 u − (δe)2 .

303

(δB)2 . Then in Eq. (A.11), two terms cancel leaving    δT : δB dv + Tr[δB]T : δB dv − δ b · δBx dm  + δ as · δBx dm > 0. (A.13) Then, since the increment in strain is constant it can be moved outside the integrals, and the final result is in terms of the average stresses and average elastic–plastic matrix and used to get the form (using normal strains only) δεε T Gδεε  0 as in the body of the paper above.

(A.10) Then Eq. (A.7) becomes      δT : δe dv + T : grad δ 2 u dv − T : (δe)2 dv   + Tr[δe]T : δe dv − δ b · δu dm   + δ as · δu dm > δ t · δu da. (A.11) Move the variation symbols inside the body force and acceleration terms, again use Eq. (A.4), and cancel the secondorder variation terms to get    δT : δe dv − T : (δe)2 dv + Tr[δe]T : δe dv   − δb · δu dm + δas · δu dm   > δt · δu da + t · δuδ(da). (A.12) Finally, to complete the loop, one can make use of the relation between the Cauchy and Piola stresses, T = J −1 TR FT and, after some lengthy algebra, return to the form of Eq. (A.1). In the body of the paper above, the displacement increment has the form δu = δBx with constant δB, so that δ 2 u = δBδx = δBδBδx. Those give that δe = δB and grad(δ 2 u) =

References Chandrasekhar, S., 1969. Ellipsoidal Figures of Equilibrium. Dover, New York. Chen, W.F., Han, D.J., 1988. Plasticity for Structural Engineers. SpringerVerlag, Berlin, New York. Dobrovolskis, A.R., 1982. Internal stresses in Phobos and other triaxial bodies. Icarus 52, 136–148. Fasso, F., Lewis, D., 2001. Stability properties of the Riemann ellipsoids. Arch. Rational Mech. Anal. 158, 259–292. Harris, A.W., 1996. The rotation rates of very small asteroids: evidence for “rubble-pile” structure. In: Proc. Lunar Planet. Sci. Conf. 27th, pp. 493– 494. Hill, R., 1958. A general theory of uniqueness and stability in elastic–plastic solids. J. Mech. Phys. Solids 6, 236–249. Hill, R., 1978. Aspects of invariance in solid mechanics. Adv. Appl. Mech. 18, 1–75. Holsapple, K.A., 2001. Equilibrium configurations of solid ellipsoidal cohesionless bodies. Icarus 154, 432–448. Lambe, T.W., Whitman, R.V., 1969. Soil Mechanics. Wiley, New York. Lebovitz, N.R., 1998. The mathematical development of the classical ellipsoids. Int. J. Engr. Sci. 36, 1407–1420. Lubarda, Vlado A., 2002. Elastoplasticity Theory. CRC Press, New York. Meirovitch, L., 1997. Principles and Techniques of Vibrations. Prentice Hall, Englewood Cliffs, NJ. Truesdell, C., Noll, W., 1965. The nonlinear field theories of mechanics. In: Flügge, S. (Ed.), Handbuch der Physik, vol. III/3. Springer-Verlag, New York. Truesdell, C., Toupin, R., 1960. The nonlinear field theories of mechanics. In: Flügge, S. (Ed.), Handbuch der Physik, vol. III/1. Springer-Verlag, New York.