Problem

If you had a long bar floating in space, what would be the compressive force at the centre of the bar, due to the self-weight of both ends?

Diagram - what is the force at point X in the middle of the bar?:

<----------------------L--------------------->, total mass M
=======================X======================   <- the bar
F---> X <---F


Summary

You should be able to simplify by cutting the bar into pieces, but that gives a different answer depending on how many pieces you use (see below). So the simplification must be wrong - but why?

My approach

Split bar in two

So, one approximation would be to cut the bar in half - two pieces of length L/2, mass M/2:

         (M/2)<-------L/2------->(M/2)
#1          X           #2   <- bar approximated as blobs #1 and #2


Force at X is G(M1.M2)/(R^2) = G (M/2)^2 / (L/2)^2 = G M^2 / L^2

Or Fx / (G. M^2 / L^2) = 1

But is that really valid? If so, shouldn't you get the same answer if you split the bar into four pieces?

Split bar into four

   (M/4)<-L/4->(M/4)<-L/4->(M/4)<-L/4->(M/4)
#1          #2    X     #3          #4


My assumption is that the force at X is the sum of the attractions of each blob on the left to every blob on the right.

Force at X = #1<>#3 + #1<>#4 + #2<>#3 + #2<>#4

('<>' being force between blobs #x and #y).

Fx / (G.M^2 / L^2) = (2/4)^-2 + (3/4)^-2 + (1/4)^-2 + (2/4)^-2 = 1.61

This is bigger than the previous result (1.61 vs 1).

Split bar into six

Similarly, if you split into 6 blobs, the total force comes out as:

Fx / (G.M^2 / L^2) = (3/6)^-2+(4/6)^-2+(5/6)^-2 + (2/6)^-2+(3/6)^-2+(4/6)^-2 + (1/6)^-2+(2/6)^-2+(3/6)^-2

Fx / (G.M^2 / L^2) = 2.00

So what's wrong with my approach? And what is the real answer?

So it seems the more pieces we split the bar into, the larger the result gets. There's clearly something wrong with my assumptions! - but what?

I'd be very glad if someone here could explain this. Thanks!

EDIT As Peter Shor pointed out, my calculations had some dodgy algebra and I'd calculated $$L^2/M^2$$ values rather than $$M^2/L^2$$. I've now corrected that - the value still increases as you divide into more masses.

I'll do a bit more work with more divisions and see if this leads to convergence or not.

-
just surround your equations with two $;) – Pratik Deoghare Nov 16 '10 at 14:55 When you separate the two pieces in the first example, each individual piece experiences the gravitational attraction from the other piece, but also a repulsion force originating from the material itself. If that wasn't the case, the material would collapse on itself. You forget to include that force in everyone of your calculations it seems to me. The total force in X is zero. – Raskolnikov Nov 16 '10 at 15:12 Clarification: as TheMachineCharmer correctly points out, the net force in the middle of the bar is zero. But what I'm trying to find is the pressure exerted on you at the middle of the bar, due to the weight of both halves of the bar. Or to ask another way, if you cut the bar in half and put a spring between the two halves, how much would the spring contract? – Sam Davies Nov 16 '10 at 15:16 Your mistake is in assuming that the center of gravity of the body will be the center of mass. That is not true. – Raskolnikov Nov 16 '10 at 15:36 @David Center of mass and center of gravity are the same in a uniform gravitational field. Center of mass is a point such that if you apply a force there, a rigid body will not rotate. Center of gravity is a point such that the torque about that point due to gravity is zero. If the gravitational field is not uniform, these could be different. In that case, there would be net torque and the rigid body would start rotating. – Mark Eichenlaub Nov 16 '10 at 21:09 3 Answers The reason you encountered higher and higher pressure at the center of the rod as you cut it into more pieces is that you were essentially approximating an integral, but the integral diverges ("is infinity" colloquially). When the rod has zero thickness, but still has mass, the density of the matter is infinite, and this leads to infinitely strong gravitational forces. To answer this question, we'll imagine the cylinder has some small, finite radius$R$. We want to find the force between the two halves of the cylinder. We'll let one half just sit stationary in space. It will create a gravitational potential. Then we'll grab the other half and pull it away to some distance$d$. The gravitational potential energy is a function of$d$. The force between the two halves of the cylinder is the derivative of the gravitational potential energy with respect to$d$when$d=0$. The problem described above is too hard. It is quite difficult to calculate the gravitational potential of a cylinder at an arbitrary point. The gravitational potential of a point mass is just$-Gm/r$, but for a cylinder that extends out in three dimensions, we need to replace$m$with the density$\rho$and then integrate over the mass of the entire cylinder. The expression for$r$, the distance from an arbitrary point outside the cylinder to a point inside it, is not very tractable. However, at a point on the axis of the cylinder, the gravitational potential is more accessible due to the extra symmetry. If we set up cylindrical coordinates with the axis of the cylinder along the z-axis, and then integrate over the bottom half of the cylinder, we get$V(z) = \int_{z'=0}^{-L/2}\int_{r=0}^{R}\int_{\theta=0}^{2\Pi} \frac{G\rho}{\sqrt{(z-z')^2+r^2}} r\textrm{d}\theta\textrm{d}r\textrm{d}z'$and doing the integral of$\theta$it's$V(z) = 2\pi G\rho\int_{z'=0}^{-L/2}\int_{r=0}^{R} \frac{1}{\sqrt{(z-z')^2+r^2}} r\textrm{d}r\textrm{d}z'$. This allows us to make an approximation. Although the half of the cylinder we use to calculate the potential must have finite width, we can calculate the potential energy by assuming that the other half of the cylinder is located perfectly along the axis. As long as the radius of the cylinder is very small compared to the length, this is a valid approximation. So the potential energy comes from integrating the previous expression for$V$along the$z$-axis for the length of the cylinder. We don't actually want the potential energy, but the derivative of the potential energy. So we imagine moving the top half up the cylinder up a little bit$dz$, and ask how the potential energy changes. Moving the entire top half of the cylinder up by$dz$is equivalent to taking a piece of thickness$dz$and slicing it off the bottom and moving it to the top. So we really just need to find the difference in the potential between the top and bottom of the top half of the cylinder and multiply by the mass-per-unit-length of the cylinder. The force between the two halves of the cylinder is$\frac{M}{L}[V(L/2) - V(0)]$That still leaves two integrals to evaluate.$V(L/2)$is easy, because it's far away from the half of the cylinder providing the gravitational potential (compared to$R$). That lets us approximate$V(L/2) = \frac{-GM}{L} \int_{-L/2}^0 \frac{1}{L/2-x}dx = -\frac{GM}{L}\ln 2$. The integral for$V(0)$is trickier, so I put it in Mathematica and got$V(0) = -\frac{GM}{L}\textrm{arcsinh}\left(\frac{L}{2R}\right)$. In the regime we are interested in ($R$small compared to$L$) the$\sinh(x)$is just$e^x/2$, so this simplifies to$V(0) = -\frac{GM}{L} \ln\left(\frac{L}{R}\right)$This gives a final answer for the force$F = \frac{M^2G}{L^2}\ln\left(\frac{L}{2R}\right)$- That's a great answer, thanks. I'm still surprised to see R in the final answer but glad to see it matches what you'd intuitively expect - ie for R << L small changes in R make little difference to the final answer (ln(x) is fairly flat for large x). – Sam Davies Nov 17 '10 at 10:09 [cont] My previous approach using integrals assumed R << L, hence r << z, so your expression for V(z) simplifies to mG int{ (1/(z-z'))dz} (where m = M/L = 2pi.R.rho = mass/unit length). This fails though because you end up with V(0) = [-ln(z')]{0,L/2} = ln(L/2) + infinity. Still not entirely sure why this doesn't work - I suspect because r << z is not true when z -> 0. – Sam Davies Nov 17 '10 at 10:19 I think that's basically the source of the problem, yeah. Glad you got something from the answer. – Mark Eichenlaub Nov 17 '10 at 10:34 As other people answered, the force is a vector quantity, so the net force at the center is indeed zero. But it seems what you have in mind is actually the pressure (stress) at that particular point. Firstly, regarding your question. What you are computing is $$\sum_{\substack{i\in\text{first part}\\j\in\text{second part}}} \mathbf F_{ij}$$ the total force acting on each particle on the first part due to the second part. This is not equivalent to "force at X". But anyway, but why do you get infinity? Because this is the answer. If we compute the gravitational field due to the 2nd part ($\lambda=M/L$): $$g(h) = G\lambda \int_0^{L/2} \frac{1}{(z+h)^2}dz = \frac{G\lambda L}{h(2h+L)},$$ we find that the gravitational force scales as$\frac1h$. Now suppose you split the first part into$N$blobs. The position of the$i$-th blob closest to the center will be$h_i=\frac{Li}{2N}$, its mass will be$\frac{\lambda L}{2N}$, so together the contribution of force of each blob will be$F_i = m_ig(h_i) \sim \frac{1/N}{i/N} = \frac1i$, and the total force is$\sim\sum_{i=1}^N\frac1i\sim\ln N$which diverges logarithmically to$\infty$. This sounds unphysical, but this is because the setting itself is unphysical — you can't have a mass with no "width" but finite mass, this makes the density infinite, and this causes the infinite force as the two part are right next to each other. So let's makes sense when the bar has size. So assume that bar is a cylinder of radius$\epsilon \to 0$. Also assume the bar's density is$\rho = M/\pi\epsilon^2L$. The gravity field along the z direction is then $$g(h) = G\rho \int_0^{L/2} \int_0^{2\pi} \int_0^{\epsilon} \frac{h+z}{((h+z)^2+r^2)^{3/2}} rdrd\phi dz = G\rho\pi \left(2 \sqrt{h^2+\epsilon ^2}-\sqrt{(2 h+L)^2+4 \epsilon ^2}+L\right),$$ and the force is now bounded above by a constant$G\rho\pi(2\epsilon+L-\sqrt{4\epsilon^2+L^2})$for varying$h$, and the total force by partitioning the 1st part would be$\sim \sum_{i=1}^N\frac{\text{constant}}N = \text{constant}$as you would expect. (Nevertheless, this constant will diverge with$1/\epsilon$as$\epsilon\to0$.) But the result above has nothing to do with the pressure. At precisely the center, there is no force, and the magnitude of force scales as$r$in its immediate neighborhood. Let's mark a small square membrane of size$a\times a\times dz$with negligible mass$a^2\rho dz$at position$z\to0above the center of the bar to record the pressure. The membrane will experience a net push from above. Since the membrane is small, we could assume it is a uniform force. The total force that pushes the membrane downwards is \begin{aligned} F &= a^2\rho dz g_{\text{total}}(z) \\ &= a^2\rho dz\cdot G\rho\pi \left(\sqrt{(L+2z)^2+4\epsilon^2}-\sqrt{(L-2z)^2+4\epsilon^2}\right) \\ &\sim \frac{4a^2z\rho dz\cdot G\rho L\pi }{\sqrt{L^2+4\epsilon^2}} \end{aligned} and thus the pressure is $$dP = -\frac F{a^2} = -\frac{4zG\rho^2 L\pi dz }{\sqrt{L^2+4\epsilon^2}} \sim 4zdz \frac{GM^2}{L^2\epsilon^4}$$ $$\implies P \sim -z^2 \frac{GM^2}{2L^2\epsilon^4}$$ note that thez$is still here. We see that the pressure is finite and maximal when we measure it precisely at the center (and decreases away from it). But this will diverge if we let$\epsilon\to0$. - beautiful...!!! – Pratik Deoghare Nov 16 '10 at 17:03 Thanks. Clearly the pressure does tend to infinity as the area on which it acts tends to zero. But what I have in mind is a "long, thin" bar - ie where the thickness of the bar is small compared to the length. What happens to your formula if you ignore the bar thickness (assume e << L), and try to calculate instead (stress * area) at the centre? – Sam Davies Nov 16 '10 at 17:06 @Kenny I think there's a mistake. Should depend on$M^2$, not$M$. I think you have calculated the force to mass ratio of a test particle at one end of the bar, not the force between the bars. – Mark Eichenlaub Nov 16 '10 at 17:47 @Mark: Actually the mistake is in the wording. I do mean the gravity force and pressure per unit mass, so it does depend only on M. – kennytm Nov 16 '10 at 17:49 @Kenny I was editing my comment while you commented. Anyway, what is pressure per unit mass? See earlier comment. – Mark Eichenlaub Nov 16 '10 at 17:52 For one thing, you've made a careless algebra mistake. You're computing$(L/M)^2$instead of$(M/L)^2\$.

-
You're right - thanks. I've corrected this in the original question now. –  Sam Davies Nov 19 '10 at 18:06
With the fix, it still won't converge, but it diverges only logarithmically, so it will diverge a lot slower now (and in fact, it will pretty much look like its converging if you go up to a moderately large number of pieces). I think Mark Eichenlaub's answer is dead on. –  Peter Shor Nov 20 '10 at 2:12