Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Archive for February, 2012

Potential due to cylindrical distribution.

Posted by peeterjoot on February 27, 2012

[Click here for a PDF of this post with nicer formatting and figures if the post had any (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

Motivation.

Consider a cylindrical distribution of mass (or charge) as in figure (\ref{fig:cylinderPotential:cylinderPotentialFig1}), with points in the cylinder given by \mathbf{r}' = (r', \theta ', z') coordinates, and the point of measurement of the potential measured at \mathbf{r} = (r, 0, 0).

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{cylinderPotentialFig1}
\caption{coordinates for evaluation of cylindrical potential.}
\end{figure}

Our potential, for a uniform distribution, will be proportional to

\begin{aligned}\begin{aligned}\phi(r) &= \int \frac{ dV'}{{\left\lvert{\mathbf{r} - \mathbf{r}'}\right\rvert}} \\ &=\int_0^R r' dr' \int_0^{2 \pi } d\theta' \int_{-L}^L \frac{dz' }{\sqrt{(z')^2+ {\left\lvert{r - r' e^{i \theta }}\right\rvert}^2}}.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.1)

Attempting to evaluate the integrals.

With

\begin{aligned}\int_{-L}^L \frac{1}{\sqrt{z^2+u^2}} \, dz=\log \left(\frac{L+\sqrt{L^2+u^2}}{-L+\sqrt{L^2+u^2}}\right)\end{aligned} \hspace{\stretch{1}}(2.2)

This is found to be

\begin{aligned}\phi(\mathbf{r}) = \int_0^R r' dr' \int_0^{2 \pi } d\theta'\log\left(\frac{L+\sqrt{L^2+ {\left\lvert{r - r' e^{ i \theta }}\right\rvert}^2}}{-L+\sqrt{L^2+ {\left\lvert{r - r' e^{ i \theta }}\right\rvert}^2}}\right)\end{aligned} \hspace{\stretch{1}}(2.3)

It is clear that we can’t evaluate this limit directly for L \rightarrow \infty since that gives us \infty/0 in the logarithm term. Presuming this can be evaluated, we must have to evaluate the complete set of integrals first, then take the limit. Based on the paper http://www.ifi.unicamp.br/~assis/J-Electrostatics-V63-p1115-1131(2005).pdf, it appears that this can be evaluated, however, the approach used therein uses mathematics a great deal more sophisticated than I can grasp without a lot of study.

Can we proceed blindly using computational tools to do the work? Attempting to evaluate the remaining integrals with Mathematica fails, since evaluation of both

\begin{aligned}\int_0^R r dr \log \left(a+\sqrt{a^2+r^2+b^2-2 r b \cos (\theta )}\right) \end{aligned} \hspace{\stretch{1}}(2.4)

and

\begin{aligned}\int_0^{2 \pi } d\theta\log \left(a+\sqrt{a^2+r^2+b^2-2 r b \cos (\theta )}\right) ,\end{aligned} \hspace{\stretch{1}}(2.5)

either time out, or take long enough that I aborted the attempt to let them evaluate.

Alternate evaluation order?

We can also attempt to evaluate this by integrating in different orders. We can for example do the r' coordinate integral first

\begin{aligned}\begin{aligned}\int_0^R &dr\frac{r}{\sqrt{a^2+r^2-2 r b \cos (\theta )+b^2}}  \\ &=\sqrt{a^2+b^2-2 b R \cos (\theta )+R^2} -\sqrt{a^2+b^2} \\ &+b \cos (\theta ) \log \left(\frac{\sqrt{a^2+b^2-2 b R \cos (\theta )+R^2}-b \cos (\theta )+R}{\sqrt{a^2+b^2}-b \cos (\theta )}\right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.6)

It should be noted that this returns a number of hard to comprehend ConditionalExpression terms, so care manipulating this expression may also be required.

If we try the angular integral first, we get

\begin{aligned}\int_0^{2 \pi } d\theta\frac{1}{\sqrt{a^2+r^2-2 r b \cos (\theta )+b^2}} =\frac{2 K\left(-\frac{4 b r}{a^2+(b-r)^2}\right)}{\sqrt{a^2+(b-r)^2}}+\frac{2 K\left(\frac{4 b r}{a^2+(b+r)^2}\right)}{\sqrt{a^2+(b+r)^2}}\end{aligned} \hspace{\stretch{1}}(2.7)

where

\begin{aligned}K(m) = F\left( \left.\frac{\pi }{2}\right| m \right)\end{aligned} \hspace{\stretch{1}}(2.8)

is the complete elliptic integral of the first kind. Actually evaluating this integral, especially in the limiting case, probably requires stepping back and thinking a bit (or a lot) instead of blindly trying to evaluate.

Posted in Math and Physics Learning. | Tagged: , , | Leave a Comment »

Potential for an infinitesimal width infinite plane. Take III

Posted by peeterjoot on February 24, 2012

[Click here for a PDF of this post with nicer formatting and figures if the post had any (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

Document generation experiment.

This little document was generated as an experiment using Mathematica and some post processing in latex.

The File menu save as latex produced latex that couldn’t be compiled, but mouse selected, copy-as latex worked out fairly well.

Post processing done included:

\begin{itemize}
\item Adding in latex prologue.
\item Stripping out the text boxes.
\item Adding in equation environments.
\item Latex generation for math output in inline text sections was uniformly poor.
\end{itemize}

Guts.

I’d like to attempt again to evaluate the potential for infinite plane distribution. The general form of our potential takes the form

\begin{aligned}\phi(\mathbf{x}) = G \rho \int \frac{1}{{{\left\lvert{\mathbf{x} - \mathbf{x}'}\right\rvert}}} dV'\end{aligned} \hspace{\stretch{1}}(2.1)

We want to evaluate this with cylindrical coordinates (r', \theta', z'), for a width \epsilon, and radius r, at distance z from the plane.

\begin{aligned}\phi (z, \epsilon , r)= 2 \pi  G \sigma  \frac{1}{\epsilon }\int _{r' = 0}^r\int _{z' = 0}^{\epsilon }\frac{r'}{\sqrt{\left(z-z'\right)^2+\left(r'\right)^2}}dz'dr'\end{aligned} \hspace{\stretch{1}}(2.2)

With the assumption that we will take the limits \epsilon \rightarrow 0, and r \rightarrow \infty. With r^2 = c/\epsilon, this does not converge. How about with r = c/\epsilon?

Performing the r’ integration (with r^2 = c/\epsilon) we find

\begin{aligned}\phi (z, \epsilon )= 2 \pi  G \sigma  \frac{1}{\epsilon }\int_{z' = 0}^{\epsilon } \left(\sqrt{\frac{c^2}{\epsilon ^2}+(z-z')^2}-\sqrt{(z-z')^2}\right) \, dz'\end{aligned} \hspace{\stretch{1}}(2.3)

Attempting to let \textit{Mathematica} evaluate this takes a long time. Long enough that I aborted the attempt to evaluate it.

Instead, first evaluating the z’ integral we have

\begin{aligned}\phi (z, \epsilon , r)=\frac{2 \pi  G \sigma }{\epsilon }\int _{r' = 0}^{c/\epsilon }\left(\log \left(\sqrt{\left(r'\right)^2+z^2}+z\right)-\log \left(\sqrt{\left(r'\right)^2+(z-\epsilon )^2}+z-\epsilon \right)\right)dr'\end{aligned} \hspace{\stretch{1}}(2.4)

This second integral can then be evaluated in reasonable time:

\begin{aligned}\begin{aligned}\phi (z, \epsilon )&= \frac{2 \pi  G \sigma }{\epsilon ^2} \left(c \log \left(\frac{\sqrt{\frac{c^2}{\epsilon ^2}+z^2}+z}{\sqrt{\frac{c^2}{\epsilon ^2}+(z-\epsilon )^2}+z-\epsilon }\right)+\epsilon  \left(z \log \left(\frac{(z-\epsilon ) \left(\sqrt{c^2+z^2 \epsilon ^2}+c\right)}{z}\right)+(\epsilon -z) \log \left(\sqrt{c^2+\epsilon ^2 (z-\epsilon )^2}+c\right)-\epsilon  \log (\epsilon  (z-\epsilon ))\right)\right) \\ &=2 \pi  G \sigma  \left(\frac{c}{\epsilon ^2}\log \left(\frac{\sqrt{c^2+z^2 \epsilon ^2}+z \epsilon }{\sqrt{c^2+\epsilon ^2(z-\epsilon )^2}+\epsilon (z-\epsilon )}\right)+\frac{z}{\epsilon } \log \left(\frac{(z-\epsilon ) \left(\sqrt{c^2+z^2 \epsilon ^2}+c\right)}{z\left(\sqrt{c^2+\epsilon ^2 (z-\epsilon )^2}+c\right)}\right)+ \log \left(\frac{\sqrt{c^2+\epsilon ^2 (z-\epsilon )^2}+c}{\epsilon  (z-\epsilon )}\right)\right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.5)

Does this have a limit as \epsilon \rightarrow 0? No, the last term is clearly divergent for c \neq 0.

Posted in Math and Physics Learning. | Tagged: , , , , | Leave a Comment »

Attempts at calculating potential distribution for infinite homogeneous plane.

Posted by peeterjoot on February 22, 2012

[Click here for a PDF of this post with nicer formatting and figures if the post had any (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

Motivation.

Part of a classical mechanics problem set was to look at what portions of momentum and angular momentum are conserved for various fields. Since this was also a previous midterm question, I’m expecting that some intuition was expected to be used determine the form of the Lagrangians, with not much effort on finding the precise form of the potentials. Here I try to calculate one such potential explicitly.

Oddly, I can explicitly calculate the potential for the infinite homogeneous plane if I start with the force and then calculate the potential, but if I start with the potential the integral also diverges? That doesn’t make any sense, so I’m wondering if I’ve miscalculated. I’ve tried a few different ways below, but can’t get a non-divergent result if I start from the integral definition of the potential instead of deriving the potential from the force after adding up all the directional force contributions (where I find of course that all the component but the one perpendicular to the plane cancel out).

Forces and potential for an infinite homogeneous plane.

Calculating the potential from the force.

For the plane, with z as the distance from the plane, and (r,\theta) coordinates in the plane, as illustrated in figure (\ref{fig:infiniteSheetPotentials:infiniteSheetPotentialsFig1}). The gravitational force from an element of mass on the plane is
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{infiniteSheetPotentialsFig1}
\caption{Coordinate choice for interaction with infinite plane mass distribution.}
\end{figure}

\begin{aligned}d\mathbf{F}&= -G \sigma m r dr d\theta \frac{z \hat{\mathbf{z}} - r \hat{\mathbf{r}}}{{\left\lvert{z \hat{\mathbf{z}} - r \hat{\mathbf{r}}}\right\rvert}^3 } \\ &= -G \sigma m r dr d\theta \frac{z \hat{\mathbf{z}} - r (\mathbf{e}_1 \cos\theta + \mathbf{e}_2 \sin\theta) }{(z^2 + r^2)^{3/2} } \\ \end{aligned}

Integrating for the total force on the test mass, noting that the sinusoidal terms vanish when integrated over a [0, 2 \pi] interval, we have

\begin{aligned}\mathbf{F} = - 2 \pi G \sigma m z \hat{\mathbf{z}} \int_0^\infty r dr \frac{1}{{(z^2 + r^2)^{3/2} }},\end{aligned} \hspace{\stretch{1}}(2.1)

a substitution r = z \tan \alpha gives us

\begin{aligned}\mathbf{F}&= - 2 \pi G \sigma m z \hat{\mathbf{z}} \int_0^{\pi/2} z \tan \alpha z \sec^2 \alpha d\alpha \frac{1}{{z^3 \sec^3 \alpha }} \\ &= - 2 \pi G \sigma m \hat{\mathbf{z}} \int_0^{\pi/2} \sin \alpha d\alpha \\ &= 2 \pi G \sigma m \hat{\mathbf{z}} (\cos(\pi/2) - \cos(0)),\end{aligned}

For

\begin{aligned}\mathbf{F} = -2 \pi G \sigma m \hat{\mathbf{z}}.\end{aligned} \hspace{\stretch{1}}(2.2)

For the Lagrangian problem we want the potential

\begin{aligned}-\boldsymbol{\nabla} \phi &= \mathbf{F} \\ -\hat{\mathbf{z}} \frac{\partial {\phi}}{\partial {z}} &= \\ \end{aligned}

or

\begin{aligned}\phi = 2 \pi G \sigma m z.\end{aligned} \hspace{\stretch{1}}(2.3)

This is a reasonable seeming answer. Our potential is of the form

\begin{aligned}\phi = m g z,\end{aligned} \hspace{\stretch{1}}(2.4)

with

\begin{aligned}g = 2 \pi G \sigma.\end{aligned} \hspace{\stretch{1}}(2.5)

Calculating the potential directly.

If we only want the potential, why start with the force? We ought to be able to work with potentials directly, and write

\begin{aligned}\phi(\mathbf{r}) = G m \rho \int \frac{dV'}{{\left\lvert{\mathbf{r}' - \mathbf{r}}\right\rvert}}.\end{aligned} \hspace{\stretch{1}}(2.6)

It’s a quick calculation to verify that is correct, and we find (provided \mathbf{r} \ne \mathbf{r}')

\begin{aligned}\mathbf{F} = - \boldsymbol{\nabla} \phi = - G m \rho \int \frac{dV' (\mathbf{r} - \mathbf{r}')}{{\left\lvert{\mathbf{r} - \mathbf{r}'}\right\rvert}^3}.\end{aligned} \hspace{\stretch{1}}(2.7)

To verify the signs, it’s helpful to refer to the diagram (\ref{fig:infiniteSheetPotentials:infiniteSheetPotentialsFig2}) which illustrates the position vectors. Now, suppose we calculate the potential directly
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{infiniteSheetPotentialsFig2}
\caption{Direction vectors for interaction with mass distribution.}
\end{figure}

Naive approach, with bogus result.

Our mass density as a function of Z is

\begin{aligned}\rho(r', \theta', z') = \lambda \delta(z')\end{aligned} \hspace{\stretch{1}}(2.8)

\begin{aligned}\phi(z)&= G m \sigma \int dz' \delta(z') \int d\theta' \int dr' \frac{r' }{((z-z')^2 + {r'}^2)^{1/2}} \\ &= 2 \pi G m \sigma \int_0^\infty \frac{r' dr' }{(z^2 + {r'}^2)^{1/2}} \\ &= 2 \pi G m \sigma z \int_0^\infty \frac{u du }{(1 + u^2)^{1/2}}.\end{aligned}

Again we can make a tangent substitution

\begin{aligned}u = \tan\alpha,\end{aligned} \hspace{\stretch{1}}(2.9)

for

\begin{aligned}\phi(z)&= 2 \pi G m \sigma z \int_0^{\pi/2} \frac{\tan \alpha \sec^2 \alpha d \alpha }{\sec\alpha} \\ &= 2 \pi G m \sigma z \int_0^{\pi/2} \tan \alpha \sec \alpha d \alpha \\ &= 2 \pi G m \sigma z \int_0^{\pi/2} \frac{\sin \alpha}{\cos^2 \alpha} d \alpha.\end{aligned}

This has the same functional form \phi = m g z as 2.4, except with

\begin{aligned}g = 2 \pi G \sigma \int_0^{\pi/2} d\alpha \frac{\sin \alpha}{\cos^2 \alpha}.\end{aligned} \hspace{\stretch{1}}(2.10)

There’s one significant and irritating difference. The integral above doesn’t have a unit value, but instead diverges (with \cos\alpha \rightarrow \infty as \alpha \rightarrow \pi/2).

What went wrong? Trouble can be seen right from the beginning. Consider the differential form above for u \gg 1

\begin{aligned}\frac{u du}{\sqrt{1 + u^2}} \approx du.\end{aligned} \hspace{\stretch{1}}(2.11)

This we are integrating on u \in [0, \infty], so long before we make the trig substitutions we are in trouble.

As the limit of a finite volume.

My guess is that we have to tie the limits of the width of the plane and its diameter, decreasing the thickness in proportion to the increase in the radius. Let’s try that.

With a finite cylinder of height \epsilon, radius R, with a measurement of the potential directly above the cylinder at height z, we have

\begin{aligned}\phi(z) = \rho G \int_0^{2\pi} d\theta' \int_0^\epsilon dz' \int_0^R r' dr' \frac{1}{{\sqrt{(z-z')^2 + {r'}^2}}}.\end{aligned} \hspace{\stretch{1}}(2.12)

Performing the \theta' integration and substituting

\begin{aligned}r' = (z - z') \tan\alpha,\end{aligned} \hspace{\stretch{1}}(2.13)

we have

\begin{aligned}\phi(z) &= 2 \pi \rho G \int_0^\epsilon dz' \int_0^{\arctan(R/(z-z'))} (z - z') \tan \alpha \sec^2 \alpha d\alpha \frac{1}{{ \sec\alpha}} \\ &= 2 \pi \rho G \int_0^\epsilon dz' \int_0^{\arctan(R/(z-z'))} (z - z') \frac{\sin \alpha}{\cos^2\alpha} d\alpha \\ &= 2 \pi \rho G \int_0^\epsilon dz' (z - z')\int_0^{\arctan(R/(z-z'))} \frac{-d\cos \alpha}{\cos^2\alpha} \\ &= 2 \pi \rho G \int_0^\epsilon dz' \evalrange{\frac{1}{{\cos\alpha}}}{0}{\arctan(R/(z-z'))} \\ &= 2 \pi \rho G \int_0^\epsilon dz' \left( \frac{1}{{\cos\arctan(R/(z-z'))}} -1 \right) \\ &= 2 \pi \rho G \int_0^\epsilon dz' \left( \sqrt{1 + \left(\frac{R}{z-z'}\right)^2} -1 \right) \\ &=2 \pi \rho G \left( -\epsilon +\int_0^\epsilon dz' \sqrt{1 + \left(\frac{R}{z-z'}\right)^2} \right)  \\ &=2 \pi \rho G \left( -\epsilon + R\int_{(z - \epsilon)/R}^{z/R} dx\sqrt{1 + \frac{1}{{x^2}}}\right)\end{aligned}

For this integral we find

\begin{aligned}\int dx \sqrt{1 + \frac{1}{{x^2}}} =\frac{x}{{\left\lvert{x}\right\rvert}}\left(\sqrt{1+x^2}+\ln(x)-\ln\left(1+\sqrt{1+x^2}\right)\right)\end{aligned} \hspace{\stretch{1}}(2.14)

Taking limits \epsilon \rightarrow 0 and R \rightarrow \infty the terms
\sqrt{1 + x^2} at either x = z/R or x = (z-\epsilon)/R tend to 1. This leaves us only with the \ln(x) contribution, so

\begin{aligned}\phi(z) =2 \pi \rho G \left( -\epsilon + R \left( \ln\left(\frac{z}{R}\right) -\ln\left(\frac{z - \epsilon}{R}\right) \right)\right)\end{aligned} \hspace{\stretch{1}}(2.15)

With \epsilon \rightarrow 0 we have the structure of a differential above, but instead of expressing the derivative in the usual forward difference form

\begin{aligned}\frac{df}{dx} = \frac{f(x + \epsilon) - f(x)}{\epsilon},\end{aligned} \hspace{\stretch{1}}(2.16)

we have to use a backwards difference, which is equivalent, provided the function f(x) is continuous at x

\begin{aligned}\frac{df}{dx} = \frac{f(x) - f(x -\epsilon)}{\epsilon}.\end{aligned} \hspace{\stretch{1}}(2.17)

We can then form the differential

\begin{aligned}df(x) = f(x) - f(x -\epsilon) = \epsilon \frac{df}{dx}.\end{aligned} \hspace{\stretch{1}}(2.18)

We also want to express the charge density \rho in terms of surface charge density \sigma, and note that these are related by \rho \Delta A \epsilon = \sigma \Delta A.

\begin{aligned}\phi(z) &=2 \pi (\rho \epsilon) G \left( -1 + R \frac{d}{dz} \ln(z/R) \right) \\ &=2 \pi \sigma G \left( -1 + \frac{d}{dz/R} \ln(z/R) \right),\end{aligned}

or

\begin{aligned}\phi(z)=2 \pi \sigma G \left( -1 + \frac{R}{z} \right)\end{aligned} \hspace{\stretch{1}}(2.19)

Looking at this result, we have the same divergent integration result as in the first attempt, and the reason for this is clear after some reflection. The limiting process for the radius and the thickness of the slice were allowed to complete independently. Before taking limits we had

\begin{aligned}\phi(z) = 2 \pi \sigma G \left( -1\frac{1}{{\epsilon}} \int_0^\epsilon dz' \sqrt{1 + \left(\frac{R}{z-z'}\right)^2} \right).\end{aligned} \hspace{\stretch{1}}(2.20)

Consider a similar, but slightly more general case, where we evaluate the limit

\begin{aligned}L = \lim_{\epsilon \rightarrow 0} \frac{1}{{\epsilon}} \int_a^{a + \epsilon} f(x) dx,\end{aligned} \hspace{\stretch{1}}(2.21)

where F'(x) = f(x), so that

\begin{aligned}L &= \lim_{\epsilon \rightarrow 0} \frac{ F(a + \epsilon) - F(a)}{\epsilon} \\ &= F'(a) \\ &= f(a).\end{aligned}

So even without evaluating the integral we expect that we’ll have

\begin{aligned}\phi(z) &= 2 \pi \sigma G \left( -1 + {\left.{{\sqrt{1 + \left(\frac{R}{z-z'}\right)^2} }}\right\vert}_{{z' = 0}} \right) \\ &\rightarrow 2 \pi \sigma G \left( -1 + \frac{R}{z} \right).\end{aligned}

This matches what was obtained in 2.19 by brute forcing the integral with Mathematica, and then evaluating the limit the hard way. Darn.

As the limit of a finite volume. Take II.

Let’s try once more. We’ll consider a homogeneous cylindrical volume of radius R, thickness \epsilon with total mass

\begin{aligned}M = \rho \pi R^2 \epsilon = \sigma \pi R^2,\end{aligned} \hspace{\stretch{1}}(2.22)

so that the area density is

\begin{aligned}\sigma = \rho \epsilon.\end{aligned} \hspace{\stretch{1}}(2.23)

Now we’ll reduce the thickness of the volume, keeping the total mass fixed, so that

\begin{aligned}\pi R^2 \epsilon = \text{constant} = \pi c^2,\end{aligned} \hspace{\stretch{1}}(2.24)

or

\begin{aligned}R = \frac{c}{\sqrt{\epsilon}}.\end{aligned} \hspace{\stretch{1}}(2.25)

We wish to evaluate

\begin{aligned}\phi(z) = 2 \pi \sigma G \frac{1}{{\epsilon}} \int_0^\epsilon dz' \int_0^{c/\sqrt{\epsilon}} r' dr' \frac{1}{{\sqrt{(z-z')^2 + {r'}^2}}}.\end{aligned} \hspace{\stretch{1}}(2.26)

Performing the integrals, (first r', then z') we find

\begin{aligned}\begin{aligned}\phi(z) = 2 \pi \sigma G \frac{1}{{2 \epsilon^2}}\Biggl(&z \left(-2 \epsilon ^2+\sqrt{\epsilon  \left(c^2+z^2 \epsilon \right)}-\sqrt{\epsilon  \left(c^2+(z-\epsilon )^2 \epsilon \right)}\right) \\ &+\epsilon  \left(\epsilon ^2+\sqrt{\epsilon  \left(c^2+(z-\epsilon )^2 \epsilon \right)}\right) \\ &+c^2 \left(-\ln \left(-z \epsilon +\sqrt{\epsilon  \left(c^2+z^2 \epsilon \right)}\right)+\ln \left(-(z - \epsilon) \epsilon +\sqrt{\epsilon  \left(c^2+(z-\epsilon )^2 \epsilon \right)}\right)\right)\Biggr)\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.27)

The \epsilon^3/\epsilon^2 term is clearly killed in the limit, and we have a -z contribution from the first term. For \epsilon \ne 0 the difference of logarithms above can be written as

\begin{aligned}\frac{-z+\epsilon +\sqrt{\frac{c^2+(z-\epsilon )^2 \epsilon }{\epsilon }}}{-z+\sqrt{\frac{c^2+z^2 \epsilon }{\epsilon }}}\end{aligned} \hspace{\stretch{1}}(2.28)

Suppose we could also validly argue that this tends to \ln(1) = 0, and the difference of square roots could also be canceled. Then we would be left with just

\begin{aligned}\phi(z) = 2 \pi \sigma G \left(-z + \frac{1}{{2\epsilon}} \sqrt{\epsilon(c^2 + (z - \epsilon)^2\epsilon} \right)\end{aligned} \hspace{\stretch{1}}(2.29)

If we also demand that z^2 \epsilon \gg c^2, then we have a z/2 contribution from the remaining square root and are left with

\begin{aligned}\phi(z) = 2 \pi \sigma G (-z/2).\end{aligned} \hspace{\stretch{1}}(2.30)

This differs from the expected result by a factor of -2, and we’ve had to do some very fishy root taking to even get that far. Employing l’H\^opital’s rule (or letting Mathematica attempt to evaluate the limit), we get infinities for the difference of logarithm term.

So, with a lot of cheating we get a result that’s similar to the expected, but not actually a match, and even to get that we had to take the limits in an invalid way. It looks like it’s back to the drawing board, but I’m not sure how to approach it.

After thinking about it a bit, perhaps the limiting process for the width needs to be explicitly accounted for using a delta function? Perhaps a QM like treatment where we express the integral in terms of some basis and look for the resolution of identity in that resolution?

Posted in Math and Physics Learning. | Tagged: , , , | Leave a Comment »

Where to find your db2dump directory for a DB2 Windows instance?

Posted by peeterjoot on February 22, 2012

On UNIX there’s a one to one mapping between instance owner userid and the path to the diaglog, and unless you’ve redirected it, you’ll find it in ~/sqllib/db2dump/. Here’s a reminder for myself of the magic required to find your db2diag.log file in a Windows instance:

E:\> db2set db2instprof
C:\ProgramData\IBM\DB2\db2build

E:\> cd /d C:\ProgramData\IBM\DB2\db2build\

C:\ProgramData\IBM\DB2\db2build> db2ilist DB2

C:\ProgramData\IBM\DB2\db2build> cd db2

C:\ProgramData\IBM\DB2\db2build\DB2> dir *log
 Volume in drive C has no label.
 Volume Serial Number is 8E04-1BBB

 Directory of C:\ProgramData\IBM\DB2\db2build\DB2

02/22/2012  03:54 PM            10,555 db2diag.log
02/22/2012  03:54 PM             1,424 db2resync.log

Posted in SQL | Tagged: , , , , | Leave a Comment »

corrupting the stack by assuming a file ulimit!

Posted by peeterjoot on February 21, 2012

I found some test code that did the following:

#if defined _AIX
   #define HANDLE_ARRAY_SIZE OPEN_MAX
#else
   #define HANDLE_ARRAY_SIZE 32768
#endif

int i = 0 ;
int fileHandleArray[HANDLE_ARRAY_SIZE] ;

memset( fileHandleArray, 0, sizeof(fileHandleArray) ) ;

do {
   fileHandleArray[i] = open( "/dev/null", O_RDWR, 0 ) ;
} while ( fileHandleArray[i++] > 0  ) ;

Fun for all!

I think this test case must have been written assuming that the ulimit for number of file handles was some nice small number. This does a nice number on the stack, and was fairly non-obvious when looking in the debugger session.

The debugger shows me the test function (testgroup2) driving a trap in a called function:

(gdb) where
#0  0x00002aaab10419d6 in gtfPrint (gtfCBInfo=0x803700008036, piFile=0x2aaab1046c9c "sqo_fileapi.C", iLine=664,
    piFormat=0x2aaab104835c "File handle value is is %d.\n") at gtfmod.h:226
#1  0x00002aaab1044035 in testgroup2 (gtfCBInfo=0x803700008036, poResult=0x803900008038) at sqo_fileapi.C:664

where my trap is on a dereference to a variable gtfCBInfo that was good earlier in the code.

        if ( gtfCBInfo->logptr == NULL )

Yet, it was clear that there were no intervening assignments to gtfCBInfo.  I was curious exactly how this corruption propagated once found, and that’s clear enough once you look at the assembly for a function call immediately after the corruption.   See all the code that sets up the parameters for the call by fetching stack spills (the loads from (%rbp)):

0x00002aaab1043fed <+1231>: mov %eax,-0xb0(%rbp)
0x00002aaab1043ff3 <+1237>: mov -0x48(%rbp),%rax
0x00002aaab1043ff7 <+1241>: lea 0x2c9e(%rip),%rdx # 0x2aaab1046c9c
0x00002aaab1043ffe <+1248>: mov $0x298,%ecx
0x00002aaab1044003 <+1253>: lea 0x4352(%rip),%rbx # 0x2aaab104835c
0x00002aaab104400a <+1260>: mov -0xb0(%rbp),%esi
0x00002aaab1044010 <+1266>: mov %rax,%rdi
0x00002aaab1044013 <+1269>: mov %esi,-0xf8(%rbp)
0x00002aaab1044019 <+1275>: mov %rdx,%rsi
0x00002aaab104401c <+1278>: mov %rcx,%rdx
0x00002aaab104401f <+1281>: mov %rbx,%rcx
0x00002aaab1044022 <+1284>: mov -0xf8(%rbp),%eax
0x00002aaab1044028 <+1290>: mov %eax,%r8d
0x00002aaab104402b <+1293>: mov $0x0,%eax
0x00002aaab1044030 <+1298>: callq 0x2aaab1041120 <_Z8gtfPrintP5gtfCBPKcmS2_z@plt>

 

Because we’ve walked all over the neighbourhood of 0(%rbp ) with the open() calls, once we try to fetch them we start passing garbage.  It’s so common to see stack corruptions manifest in even nastier ways (like a trap on return from a function call because we’ve wiped out the return address), so I was a bit suprised to see such a run of the mill seeming corruption here as a result.  Perhaps if we hadn’t called any functions, we would have trapped on unwind instead?

Posted in C/C++ development and debugging. | Tagged: , | Leave a Comment »

Compilation of class notes for phy454h1s, continuum mechanics (so far).

Posted by peeterjoot on February 17, 2012

Have collected all my pre-midterm continuum mechanics notes into a single document. The individual pdfs below are still available, but won’t be updated further.

Feb 17, 2012 Flow in a pipe. Gravity driven flow of a film.

Feb 15, 2012 Worked examples of unidirectional Navier-Stokes solutions.

Feb 10, 2012 Navier-Stokes equation.

Feb 8, 2012 Newtonian fluids. Mass conservation. Constitutive relation. Incompressible fluids.

Feb 8, 2012 Problem Set 1. Stress, Strain, Traction vector. Force free equilibrium.

Feb 3, 2012 Phasor description of elastic waves. Fluid dynamics.

Feb 1, 2012 P-waves and S-waves.

Jan 27, 2012 Compatibility condition and elastostatics.

Jan 25, 2012 Constitutive relationship.

Jan 20, 2012 Strain tensor components.

Jan 18, 2012 Strain tensor review. Stress tensor.

Jan 13, 2012 Introduction and strain tensor.

Jan 11, 2012 Overview.

Posted in Math and Physics Learning. | Tagged: , , , , , , , , , , , , | 1 Comment »

PHY454H1S Continuum Mechanics. Lecture 12: Flow in a pipe. Gravity driven flow of a film. Taught by Prof. K. Das.

Posted by peeterjoot on February 17, 2012

[Click here for a PDF of this post with nicer formatting and figures if the post had any (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

Disclaimer.

Peeter’s lecture notes from class. May not be entirely coherent.

Review. Steady rectilinear flow.

Steady:

\begin{aligned}\frac{\partial {}}{\partial {t}} = 0\end{aligned} \hspace{\stretch{1}}(2.1)

Rectilinear is a unidirectional flow such as

\begin{aligned}\mathbf{u} = \hat{\mathbf{x}} u( x, y, z ),\end{aligned} \hspace{\stretch{1}}(2.2)

\begin{enumerate}
\item
Utilizing an incompressibility assumption \boldsymbol{\nabla} \cdot \mathbf{u} = 0, so for this case we have

\begin{aligned}\frac{\partial {u}}{\partial {x}} = 0\end{aligned}

or

\begin{aligned}u = u(y, z)\end{aligned}

Note that Prof. Das called this a continuity requirement, and justified this label with the relation

\begin{aligned}\frac{d\rho}{dt} = \rho (\boldsymbol{\nabla} \cdot \mathbf{u}),\end{aligned} \hspace{\stretch{1}}(2.3)

which was a consequence of mass conservation. It’s still not clear to me why he would call this a continuity requirement.

\item Nonlinear term is zero. (\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u} = 0
\item p = p(x). Since \frac{d^2 p}{dx^2} = 0 we also have \frac{dp}{dx} = -G, a constant.

\item \mu \left( \frac{\partial^2 {{u}}}{\partial {{y}}^2} + \frac{\partial^2 {{u}}}{\partial {{z}}^2} \right) = G

\end{enumerate}

Solution by intuition.

Two examples that we have solved analytically are illustrated in figure (\ref{fig:continuumL12:continuumL12fig1}) and figure (\ref{fig:continuumL12:continuumL12fig2})

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig1}
\caption{Simple shear flow}
\end{figure}
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig2}
\caption{Channel flow}
\end{figure}

Sometimes we can utilize solutions already found to understand the behaviour of more complex systems. Combining the two we can look at flow over a plate as in figure (\ref{fig:continuumL12:continuumL12fig3})

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig3}
\caption{Flow on a plate}
\end{figure}

Example 2. Fluid in a container. If the surface tension is altered on one side, we induce a flow on the surface, leading to a circulation flow. This can be done for example, by introducing a heat source or addition of surfactant.

This is illustrated in figure (\ref{fig:continuumL12:continuumL12fig4})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig4}
\caption{Circulation flow induced by altering surface tension.}
\end{figure}

This sort of flow is hard to analyze, only first done by Steve Davis in the 1980’s. The point here is that we can use some level of intuition to guide our attempts at solution.

Flow down a pipe.

Reading: section 2 from [1].

Recall that the Navier-Stokes equation is

\begin{aligned}\boxed{\rho \frac{\partial {\mathbf{u}}}{\partial {t}} + \rho (\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u} = - \boldsymbol{\nabla} p + \mu \boldsymbol{\nabla}^2 \mathbf{u} + \rho \mathbf{f}.}\end{aligned} \hspace{\stretch{1}}(4.4)

We need to express this in cylindrical coordinates (r, \theta, z) as in figure (\ref{fig:continuumL12:continuumL12fig5})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig5}
\caption{Flow through a pipe.}
\end{figure}

Our gradient is

\begin{aligned}\boldsymbol{\nabla} = \hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \frac{\hat{\boldsymbol{\theta}}}{r} \frac{\partial {}}{\partial {\theta}} + \hat{\mathbf{z}} \frac{\partial {}}{\partial {z}},\end{aligned} \hspace{\stretch{1}}(4.5)

For our Laplacian we find

\begin{aligned}\boldsymbol{\nabla}^2 &= \left(\hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \frac{\hat{\boldsymbol{\theta}}}{r} \frac{\partial {}}{\partial {\theta}} + \hat{\mathbf{z}} \frac{\partial {}}{\partial {z}}\right) \cdot\left(\hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \frac{\hat{\boldsymbol{\theta}}}{r} \frac{\partial {}}{\partial {\theta}} + \hat{\mathbf{z}} \frac{\partial {}}{\partial {z}}\right) \\ &=\partial_{rr} + \frac{\hat{\boldsymbol{\theta}}}{r} \cdot (\partial_\theta \hat{\mathbf{r}}) \partial_r+ \frac{1}{{r}} \partial_\theta \left( \frac{1}{{r}} \partial_\theta \right)+ \partial_{zz} \\ &=\partial_{rr} + \frac{1}{{r}} \partial_r + \frac{1}{{r^2}} \partial_{\theta\theta} + \partial_{zz},\end{aligned}

which we can write as

\begin{aligned}\boldsymbol{\nabla}^2 = \frac{1}{{r}} \frac{\partial {}}{\partial {r}} \left( r \frac{\partial {}}{\partial {r}} \right) + \frac{1}{{r^2}} \frac{\partial^2 {{}}}{\partial {{\theta}}^2} + \frac{\partial^2 {{}}}{\partial {{z}}^2}.\end{aligned} \hspace{\stretch{1}}(4.6)

NS takes the form

\begin{aligned}\boxed{\begin{aligned}\rho \frac{\partial {\mathbf{u}}}{\partial {t}} &+ \rho \left(u_r \frac{\partial {}}{\partial {r}} + \frac{u_\theta}{r} \frac{\partial {}}{\partial {\theta}} + u_z \frac{\partial {}}{\partial {z}} \right) \mathbf{u} =  \\ &- \left(\hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \frac{\hat{\boldsymbol{\theta}}}{r} \frac{\partial {}}{\partial {\theta}} + \hat{\mathbf{z}} \frac{\partial {}}{\partial {z}}\right)p + \mu \left(\frac{1}{{r}} \frac{\partial {}}{\partial {r}} \left( r \frac{\partial {}}{\partial {r}} \right) + \frac{1}{{r^2}} \frac{\partial^2 {{}}}{\partial {{\theta}}^2} + \frac{\partial^2 {{}}}{\partial {{z}}^2} \right)\mathbf{u} + \rho \mathbf{f}.\end{aligned}}\end{aligned} \hspace{\stretch{1}}(4.7)

For steady state and incompressible fluids in the absence of body forces we have

\begin{aligned}\left(\hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \frac{\hat{\boldsymbol{\theta}}}{r} \frac{\partial {}}{\partial {\theta}} + \hat{\mathbf{z}} \frac{\partial {}}{\partial {z}}\right)p = \mu \left(\frac{1}{{r}} \frac{\partial {}}{\partial {r}} \left( r \frac{\partial {}}{\partial {r}} \right) + \frac{1}{{r^2}} \frac{\partial^2 {{}}}{\partial {{\theta}}^2} + \frac{\partial^2 {{}}}{\partial {{z}}^2} \right)\mathbf{u},\end{aligned} \hspace{\stretch{1}}(4.8)

or, in coordinates

\begin{aligned}\frac{\partial {p}}{\partial {r}}  &= \mu \left(\frac{1}{{r}} \frac{\partial {}}{\partial {r}} \left( r \frac{\partial {}}{\partial {r}} \right) + \frac{1}{{r^2}} \frac{\partial^2 {{}}}{\partial {{\theta}}^2} + \frac{\partial^2 {{}}}{\partial {{z}}^2} \right)u_r \\ \frac{1}{r} \frac{\partial {p}}{\partial {\theta}}&= \mu \left(\frac{1}{{r}} \frac{\partial {}}{\partial {r}} \left( r \frac{\partial {}}{\partial {r}} \right) + \frac{1}{{r^2}} \frac{\partial^2 {{}}}{\partial {{\theta}}^2} + \frac{\partial^2 {{}}}{\partial {{z}}^2} \right)u_\theta \\ \frac{\partial {p}}{\partial {z}}&= \mu \left(\frac{1}{{r}} \frac{\partial {}}{\partial {r}} \left( r \frac{\partial {}}{\partial {r}} \right) + \frac{1}{{r^2}} \frac{\partial^2 {{}}}{\partial {{\theta}}^2} + \frac{\partial^2 {{}}}{\partial {{z}}^2} \right)u_z\end{aligned} \hspace{\stretch{1}}(4.9)

With an assumption that we have no radial or circulatory flows (u_r = u_\theta = 0), and with u_z = w assumed to only have a radial dependence, our velocity is

\begin{aligned}\mathbf{u} = \hat{\mathbf{z}} w(r),\end{aligned} \hspace{\stretch{1}}(4.12)

and an assumption of linear pressure dependence

\begin{aligned}\frac{dp}{dz} = -G,\end{aligned} \hspace{\stretch{1}}(4.13)

then NS takes the final simple form

\begin{aligned}\frac{1}{{r}} \frac{d}{dr} \left( r \frac{dw}{dr} \right) = - \frac{G}{\mu}.\end{aligned} \hspace{\stretch{1}}(4.14)

Solving this we have

\begin{aligned}r \frac{dw}{dr} = - \frac{G r^2}{2\mu} + A\end{aligned} \hspace{\stretch{1}}(4.15)

\begin{aligned}w = -\frac{G r^2}{4 \mu} + A \ln(r) + B\end{aligned} \hspace{\stretch{1}}(4.16)

Requiring finite solutions for r = 0 means that we must have A = 0. Also w(a) = 0, we have B = G a^2/4 \mu so we must have

\begin{aligned}w(r) = \frac{G}{4 \mu}( a^2 - r^2 )\end{aligned} \hspace{\stretch{1}}(4.17)

Example: Gravity driven flow of a liquid film

(This is one of our Professor’s favorite problems).

Coordinates as in figure (\ref{fig:continuumL12:continuumL12fig6})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig6}
\caption{Gravity driven flow down an inclined plane.}
\end{figure}

\begin{aligned}\mathbf{u} = \hat{\mathbf{x}} u(y)\end{aligned} \hspace{\stretch{1}}(5.18)

Boundary conditions

\begin{enumerate}
\item u(y = 0) = 0
\item Tangential stress at the air-liquid interface y = h is equal.

\begin{aligned}\boldsymbol{\tau} \cdot (\boldsymbol{\sigma}_l \cdot \hat{\mathbf{n}}) = \boldsymbol{\tau} \cdot (\boldsymbol{\sigma}_a \cdot \hat{\mathbf{n}}),\end{aligned} \hspace{\stretch{1}}(5.19)

\end{enumerate}

We write

\begin{aligned}\boldsymbol{\tau} &= \begin{bmatrix}1 \\ 0 \\ 0\end{bmatrix} \\ \hat{\mathbf{n}} &= \begin{bmatrix}0 \\ 1 \\ 0\end{bmatrix} \\ \end{aligned} \hspace{\stretch{1}}(5.20)

and seek simultaneous solutions to the pair of stress tensor equations

\begin{aligned}\sigma_{ij}^l &= - p \delta_{ij} + \mu^l \left( \frac{\partial {u_i}}{\partial {x_j}} +\frac{\partial {u_j}}{\partial {x_i}}\right) \\ \sigma_{ij}^a &= - p \delta_{ij} + \mu^a \left( \frac{\partial {u_i}}{\partial {x_j}} +\frac{\partial {u_j}}{\partial {x_i}}\right).\end{aligned} \hspace{\stretch{1}}(5.23)

In general this requires an iterated approach, solving for one with an initial approximation of the other, then switching and tuning the numerical method carefully for convergence.

We expect that the flow of liquid will induce a flow of air at the interface, but may be able to make a one-sided approximation. Let’s see how far we get before we have to introduce any approximations and compute the traction vector for the liquid

\begin{aligned}\boldsymbol{\sigma}^l \cdot \hat{\mathbf{n}} &= \begin{bmatrix}-p & \mu^l {\partial {u}}/{\partial {y}} & 0 \\ \mu^l {\partial {u}}/{\partial {y}} & -p & 0 \\ 0 & 0 & 0\end{bmatrix}\begin{bmatrix}0 \\ 1 \\ 0\end{bmatrix} \\ &=\begin{bmatrix}\mu^l {\partial {u}}/{\partial {y}} \\ -p \\ 0\end{bmatrix}\end{aligned}

So

\begin{aligned}\boldsymbol{\tau} \cdot (\boldsymbol{\sigma}^l \cdot \hat{\mathbf{n}})=\begin{bmatrix}1 & 0 & 0\end{bmatrix}\begin{bmatrix}\mu^l {\partial {u}}/{\partial {y}} \\ -p \\ 0\end{bmatrix}=\mu^l \frac{\partial {u}}{\partial {y}}\end{aligned} \hspace{\stretch{1}}(5.25)

Our boundary value condition is therefore

\begin{aligned}{\left.{{\mu^l \frac{\partial {u^l}}{\partial {y}}}}\right\vert}_{{y = h}} ={\left.{{\mu^a \frac{\partial {u^a}}{\partial {y}}}}\right\vert}_{{y = h}}\end{aligned} \hspace{\stretch{1}}(5.26)

When can we decouple this, treating only the liquid? Observe that we have

\begin{aligned}{\left.{{\frac{\partial {u^l}}{\partial {y}}}}\right\vert}_{{y = h}} ={\left.{{\frac{\mu^a}{\mu^l} \frac{\partial {u^a}}{\partial {y}}}}\right\vert}_{{y = h}}\end{aligned} \hspace{\stretch{1}}(5.27)

so if

\begin{aligned}\frac{\mu_a}{\mu_l} \ll 1\end{aligned} \hspace{\stretch{1}}(5.28)

we can treat only the liquid portion of the problem, with a boundary value condition

\begin{aligned}{\left.{{\frac{\partial {u^l}}{\partial {y}}}}\right\vert}_{{y = h}} = 0.\end{aligned} \hspace{\stretch{1}}(5.29)

Let’s look at the component of the traction vector in the direction of the normal (liquid pressure acting on the air)

\begin{aligned}\hat{\mathbf{n}} \cdot (\boldsymbol{\sigma}^l \cdot \hat{\mathbf{n}}) = \hat{\mathbf{n}} \cdot (\boldsymbol{\sigma}^a \cdot \hat{\mathbf{n}}) \end{aligned} \hspace{\stretch{1}}(5.30)

or

\begin{aligned}\begin{bmatrix}0 & 1 & 0\end{bmatrix}\begin{bmatrix}\mu^l \frac{\partial {u}}{\partial {y}} \\ -p^l \\ 0\end{bmatrix}= -{\left.{{p^l}}\right\vert}_{{y = h}} = -{\left.{{p^a}}\right\vert}_{{y = h}}\end{aligned} \hspace{\stretch{1}}(5.31)

i.e. We have pressure matching at the interface.

Our body force is

\begin{aligned}\mathbf{f} = \begin{bmatrix}g \sin\alpha \\ -g \cos\alpha \\ 0\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(5.32)

Referring to the Navier-Stokes equation 4.4, we see that our only surviving parts are

\begin{subequations}

\begin{aligned}0 = -\frac{\partial {p}}{\partial {x}} + \mu \frac{\partial^2 {{u}}}{\partial {{y}}^2} + \rho g \sin\alpha \end{aligned} \hspace{\stretch{1}}(5.33a)

\begin{aligned}0 = -\frac{\partial {p}}{\partial {y}} - \rho g \cos\alpha \end{aligned} \hspace{\stretch{1}}(5.33b)

\begin{aligned}0 = -\frac{\partial {p}}{\partial {z}} \end{aligned} \hspace{\stretch{1}}(5.33c)

\end{subequations}

The last gives us p \ne p(z). Integrating the second we have

\begin{aligned}p = \rho g y \cos\alpha + p_1\end{aligned} \hspace{\stretch{1}}(5.34)

Since p = p_{\text{atm}} at y = h, we have

\begin{aligned}p_{\text{atm}} = \rho g h \cos\alpha + p_1\end{aligned} \hspace{\stretch{1}}(5.35)

Our first NS equation 5.33a becomes

\begin{aligned}0 = \mu \frac{\partial^2 {{u}}}{\partial {{y}}^2} + g \sin\alpha,\end{aligned} \hspace{\stretch{1}}(5.36)

or

\begin{aligned}\frac{\partial^2 {{u}}}{\partial {{y}}^2} = -\frac{g}{\mu} \sin\alpha\end{aligned} \hspace{\stretch{1}}(5.37)

Solving we have

\begin{aligned}u = - \rho g \frac{\sin\alpha}{2 \mu} y^2 + A y + B\end{aligned} \hspace{\stretch{1}}(5.38)

With

\begin{aligned}u(0) &= 0 \\ {\left.{{\frac{\partial {u}}{\partial {y}}}}\right\vert}_{{y = h}} &= 0\end{aligned} \hspace{\stretch{1}}(5.39)

\begin{aligned}u = \rho g \frac{\sin\alpha}{2 \mu} \left( 2 h y - y^2 \right) .\end{aligned} \hspace{\stretch{1}}(5.41)

This velocity distribution is illustrated figure (\ref{fig:continuumL12:continuumL12fig7}).

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig7}
\caption{Velocity streamlines for flow down a plane.}
\end{figure}

It’s important to note that in these problems we have to derive our boundary value conditions! They are not given.

In this discussion, the height h was assumed to be constant, with the tangential direction constant and parallel to the surface that the liquid is flowing on. It’s claimed in class that this is actually a consequence of surface tension only! That’s not at all intuitive, but will be covered when we learn about “stability conditions”.

Study note.

Memorizing the NS equation is required for midterm, but more complex stuff (like cylindrical forms of the strain tensor if required) will be given.

References

[1] D.J. Acheson. Elementary fluid dynamics. Oxford University Press, USA, 1990.

Posted in Math and Physics Learning. | Tagged: , , , , , , , , , , | Leave a Comment »

PHY454H1S Continuum Mechanics. Lecture 11: Worked examples of Navier-Stokes solutions. Taught by Prof. K. Das.

Posted by peeterjoot on February 16, 2012

[Click here for a PDF of this post with nicer formatting and figures if the post had any (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

Disclaimer.

Peeter’s lecture notes from class. May not be entirely coherent.

Navier-Stokes equation.

The Navier-Stokes equation (our fluids equivalent to Newton’s second law) was found to be

\begin{aligned}\rho \frac{\partial {\mathbf{u}}}{\partial {t}} + \rho (\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u} = - \boldsymbol{\nabla} p + \mu \boldsymbol{\nabla}^2 \mathbf{u} + \rho \mathbf{f}.\end{aligned} \hspace{\stretch{1}}(2.1)

In this course we’ll focus on the incompressible case where we have

\begin{aligned}\boldsymbol{\nabla} \cdot \mathbf{u} = 0\end{aligned} \hspace{\stretch{1}}(2.2)

We watched a video of the rocking tank as in figure (\ref{fig:continuumL11:continuumL11fig1}). The boundary condition that accounted for the matching of the die marker is that we have no slipping at the interface. Writing \boldsymbol{\tau} for the tangent to the interface then this condition at the interface is described mathematically by the conditions

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL11fig1}
\caption{Rocking tank velocity matching.}
\end{figure}

\begin{aligned}\mathbf{u}_A \cdot \boldsymbol{\tau} &= \mathbf{u}_B \cdot \boldsymbol{\tau} \\ \mathbf{u}_A \cdot \hat{\mathbf{n}} &= \mathbf{u}_B \cdot \hat{\mathbf{n}}.\end{aligned} \hspace{\stretch{1}}(2.3)

Referring to figure (\ref{fig:continuumL11:continuumL11fig2}) where the tangents and normals are depicted an example representation of the normal and tangent vectors for the fluids are

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL11fig2}
\caption{Normals and tangents at interface for 2D system}
\end{figure}

\begin{aligned}\boldsymbol{\tau} &= \begin{bmatrix}1 \\ 0\end{bmatrix} \\ \hat{\mathbf{n}} &= \begin{bmatrix}0 \\ 1\end{bmatrix} \end{aligned} \hspace{\stretch{1}}(2.5)

For the traction vector

\begin{aligned}T_i = \sigma_{ij} n_j,\end{aligned} \hspace{\stretch{1}}(2.7)

we also have at the interface we must have matching of

\begin{aligned}\boldsymbol{\tau} \cdot \mathbf{T}.\end{aligned} \hspace{\stretch{1}}(2.8)

More explicitly, in coordinates this is

\begin{aligned}{\left.{{\tau_i (\sigma_{ij} n_j)}}\right\vert}_{{A}} ={\left.{{\tau_i (\sigma_{ij} n_j)}}\right\vert}_{{B}}\end{aligned} \hspace{\stretch{1}}(2.9)

Steady incompressible rectilinear (unidirectional) flow.

In this case we can fix our axis so that

\begin{aligned}\mathbf{u} = \hat{\mathbf{x}} u(x, y, z, t),\end{aligned} \hspace{\stretch{1}}(3.10)

where the velocity components in the other directions

\begin{aligned}v &= 0 \\ w &= 0\end{aligned} \hspace{\stretch{1}}(3.11)

are both zero. Symbolically, the steady state condition is

\begin{aligned}\frac{\partial {\mathbf{u}}}{\partial {t}} = 0.\end{aligned} \hspace{\stretch{1}}(3.13)

We start with the incompressibility condition, which written explicitly, is

\begin{aligned}\boldsymbol{\nabla} \cdot \mathbf{u} = 0,\end{aligned} \hspace{\stretch{1}}(3.14)

or

\begin{aligned}\frac{\partial {u}}{\partial {x}} + \frac{\partial {v}}{\partial {y}} + \frac{\partial {w}}{\partial {z}} = 0\end{aligned} \hspace{\stretch{1}}(3.15)

This implies

\begin{aligned}\frac{\partial {u}}{\partial {x}} = 0\end{aligned} \hspace{\stretch{1}}(3.16)

so our velocity can only be function of the y and z coordinates only

\begin{aligned}u = u(y, z).\end{aligned} \hspace{\stretch{1}}(3.17)

The non-linear term of the Navier-Stokes equation takes the form

\begin{aligned}(\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}&=\left(u \frac{\partial {}}{\partial {x}}+\not{{v}} \frac{\partial {}}{\partial {y}}+\not{{w}} \frac{\partial {}}{\partial {z}} \right) (\hat{\mathbf{x}} u( y, z) + \hat{\mathbf{y}} \not{{v}} + \hat{\mathbf{z}} \not{{w}} ) \\ &= \hat{\mathbf{x}} u \not{{\frac{\partial {u}}{\partial {x}}}} \\ &= 0.\end{aligned}

With incompressibility and u = v = 0 conditions killing this term, and the steady state condition 3.13 killing the \rho {\partial {\mathbf{u}}}/{\partial {t}} term, the Navier-Stokes equation for this incompressible unidirectional steady state flow (in the absence of body forces) is reduced to

\begin{aligned}0 = - \boldsymbol{\nabla} p + \mu \boldsymbol{\nabla}^2 \mathbf{u}.\end{aligned} \hspace{\stretch{1}}(3.18)

In coordinates this is

\begin{subequations}

\begin{aligned}\frac{\partial {p}}{\partial {x}} = \mu \left( \frac{\partial^2 {{u}}}{\partial {{y}}^2} + \frac{\partial^2 {{u}}}{\partial {{z}}^2} \right) \end{aligned} \hspace{\stretch{1}}(3.19a)

\begin{aligned}\frac{\partial {p}}{\partial {y}} = 0 \end{aligned} \hspace{\stretch{1}}(3.19b)

\begin{aligned}\frac{\partial {p}}{\partial {z}} = 0.\end{aligned} \hspace{\stretch{1}}(3.19c)

\end{subequations}

Operating on the first with an x partial we find

\begin{aligned}\frac{\partial^2 {{p}}}{\partial {{x}}^2} = \mu \left( \frac{\partial^2 {{}}}{\partial {{y}}^2} \not{{\frac{\partial {u}}{\partial {x}}}} + \frac{\partial^2 {{}}}{\partial {{z}}^2} \not{{ \frac{\partial {u}}{\partial {x}} }} \right) = 0\end{aligned} \hspace{\stretch{1}}(3.20)

Since we have

\begin{aligned}\frac{\partial^2 {{p}}}{\partial {{x}}^2} = 0\end{aligned} \hspace{\stretch{1}}(3.21)

we also have

\begin{aligned}\frac{d^2 p}{dx^2} = 0,\end{aligned} \hspace{\stretch{1}}(3.22)

so our pressure must be linear with position

\begin{aligned}p = A x + B,\end{aligned} \hspace{\stretch{1}}(3.23)

as illustrated in figure (\ref{fig:continuumL11:continuumL11fig3})

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL11fig3}
\caption{Pressure gradient in 1D system.}
\end{figure}

\begin{aligned}p =\left\{\begin{array}{l l}p_0 & \quad \mbox{latex x = 0$} \\ p_L & \quad \mbox{x = L} \end{array}\right.\end{aligned} \hspace{\stretch{1}}(3.24)$

we have

\begin{aligned}p = \frac{p_L - p_0}{L} x + p_0\end{aligned} \hspace{\stretch{1}}(3.25)

and

\begin{aligned}\frac{dp}{dx} = \frac{p_L - p_0}{L} = \text{constant} \equiv -G\end{aligned} \hspace{\stretch{1}}(3.26)

Example: Shearing flow.

The flows of this sort don’t have to be trivial. For example, even with constant pressure (p_0 = p_L) as in figure (\ref{fig:continuumL11:continuumL11fig4}) we can have a “shearing flow” where the fluids at the top surface are not necessarily moving at the same rates as the fluid below that surface. We have fluid flow in the x direction only, and our velocity is a function only of the y coordinate.

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL11fig4}
\caption{Velocity variation with height in shearing flow.}
\end{figure}

\begin{aligned}\mathbf{u} &= \hat{\mathbf{x}} u(y) \\ G &= 0 \\ u(0) &= 0 \\ u(h) &= U.\end{aligned} \hspace{\stretch{1}}(3.27)

For such a flow 3.19a simplifies to

\begin{aligned}\frac{d^2 u}{dy^2} = 0\end{aligned} \hspace{\stretch{1}}(3.31)

with solution

\begin{aligned}u = \frac{U}{h} y + u(0) = \frac{U}{h} y.\end{aligned} \hspace{\stretch{1}}(3.32)

Example: Channel flow

\begin{aligned}\mathbf{u} &= \hat{\mathbf{x}} u(y) \\ G &= - \frac{dp}{dx} \ne 0\end{aligned} \hspace{\stretch{1}}(3.33)

This time our simplified Navier-Stokes equation 3.19a is reduced to something slightly more complicated

\begin{aligned}\mu \frac{d^2 u}{dy^2} = -G,\end{aligned} \hspace{\stretch{1}}(3.35)

with solution

\begin{aligned}u = -\frac{G}{2 \mu} y^2 + A y + B.\end{aligned} \hspace{\stretch{1}}(3.36)

The boundary value conditions with the coordinate system in use illustrated in figure (\ref{fig:continuumL11:continuumL11fig5}) require the velocity to be zero at the interface (the pipe walls preventing flow in the interior of the pipe)

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL11fig5}
\caption{1D Channel flow coordinate system setup.}
\end{figure}

\begin{aligned}u(\pm h) = -\frac{G}{2 \mu} h^2 \pm A h + B = 0\end{aligned} \hspace{\stretch{1}}(3.37)

One solution, immediately evident is,

\begin{aligned}A &= 0 \\ B &= \frac{G}{2 \mu} h^2,\end{aligned} \hspace{\stretch{1}}(3.38)

so our solution becomes

\begin{aligned}u = \frac{G}{2 \mu} \Bigl( h^2 - y^2 \Bigr),\end{aligned} \hspace{\stretch{1}}(3.40)

a parabolic velocity flow. This is illustrated graphically in figure (\ref{fig:continuumL11:continuumL11fig6}).

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL11fig6}
\caption{Parabolic velocity distribution.}
\end{figure}

It’s clear that this is maximized by y = 0, but we can also see this by computing

\begin{aligned}\frac{du}{dy} = \frac{G}{\mu} y = 0.\end{aligned} \hspace{\stretch{1}}(3.41)

This maximum is

\begin{aligned}u_{\text{max}} = \frac{G}{2\mu} h^2\end{aligned} \hspace{\stretch{1}}(3.42)

The flux, or flow rate is

\begin{aligned}Q &= \iint_S \mathbf{u} \cdot \hat{\mathbf{x}} ds \\ &= \int_0^1 dz \int_{-h}^h dy u(y) \\ &=\frac{2 G h^3}{3}\end{aligned}

Let’s now compute the strain (e_{ij}) and the stress (\sigma_{ij} = -p \delta_{ij} + 2 \mu e_{ij})

\begin{aligned}e_{12} &= e_{21} = \frac{1}{{2}} \left( \frac{\partial {u}}{\partial {y}} \right) = - \frac{G y}{2 \mu} \\ e_{11} &= \frac{\partial {u}}{\partial {x}} = 0 \\ e_{22} &= \frac{\partial {v}}{\partial {y}} = 0\end{aligned} \hspace{\stretch{1}}(3.43)

stress

\begin{aligned}\sigma_{12} = 2 \mu e_{12} = -G y\end{aligned} \hspace{\stretch{1}}(3.46)

This can be used to compute the forces on the inner surfaces of the tube. As illustrated in figure (\ref{fig:continuumL11:continuumL11fig7}), our normals at \pm h are \mp \hat{\mathbf{y}} respectively. The traction vector in the y direction is at y = h is

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL11fig7}
\caption{Normals in 1D channel flow system}
\end{figure}

\begin{aligned}T_i = \sigma_{i 2} {\left.{{n_2}}\right\vert}_{{y = h}} = G h,\end{aligned} \hspace{\stretch{1}}(3.47)

so that

\begin{aligned}\mathbf{T} = \hat{\mathbf{x}} G h\end{aligned} \hspace{\stretch{1}}(3.48)

(here the x directionality comes from the i = 1 index of the stress tensor).

\begin{aligned}F_{x_L} = \hat{\mathbf{x}} \cdot \mathbf{T} = G h\end{aligned} \hspace{\stretch{1}}(3.49)

The total force is then

\begin{aligned}\int_0^L F_x dx = + G h L\end{aligned} \hspace{\stretch{1}}(3.50)

FIXME: ask in class. This is the tangential force at the boundary of the wall. What is it a force on? If it is tangential, how can it act on the wall? It could act on an impediment placed right up next to the wall, but if that’s the case, why are we integrating from x = 0 to x = L?

Posted in Math and Physics Learning. | Tagged: , , , , , , , , , | Leave a Comment »

A workaround for the viral effect of incorrect Lotus notes automatic add to contacts (bad distribution lists)

Posted by peeterjoot on February 14, 2012

In Lotus Notes version 8, once an email is sent with an incorrect email address, it can have an almost viral effect. For example, we have a couple internal distribution lists that have been used incorrectly a few times:

askOSS@torolab.ibm.com
askBPS@torolab.ibm.com

(The actual distribution lists have canonical lotus notes addresses askOSS/…/…@…)

I’m on one such distribution list, and this incorrect address has been repeatedly added to my local address book, so once I try to send an email to our askOSS list, it bounces.

The workaround is to search through the local address book, find the bad addresses, and delete them. This unfortunately has to be done on any replicas separately, and it’s temporary because somebody ends up not knowing to do this, and the bad address ends up back in everybodies local address books again.

I asked our IT support folks about this. They’ve opened a defect against Lotus notes for this issue, but have also provided the following workaround:

As a temporary fix

  1. You will have to delete the incorrect address from your address book under recent contacts which you have been doing
  2. create a contact in the address book with the incorrect address as the main contact name with the right mail address in the contents of it
  3. you can disable the recent contact feature to add the names to the recent contacts field (before doing this you will have clear all the bad address from Contacts|recent contacts field)

Steps to disable recent contacts:

  • Click File|preferences|contacts
  • Check “Do not automatically add contacts to the recent contacts view”
  • click ok

I’m hopeful that this will deal with the problem.  The side effect will be that I will have a smaller contact list, but I think I can live with that.

Posted in Development environment | Tagged: | Leave a Comment »

Runge-Lenz vector conservation

Posted by peeterjoot on February 13, 2012

[Click here for a PDF of this post with nicer formatting and figures if the post had any (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

Motivation.

Notes from Prof. Poppitz’s phy354 classical mechanics lecture on the Runge-Lenz vector, a less well known conserved quantity for the 3D 1/r potentials that can be used to solve the Kepler problem.

Motivation: The Kepler problem.

We can plug away at the Lagrangian in cylindrical coordinates and find eventually

\begin{aligned}\int_{\phi_0}^\phi d\phi = \int_{r_0}^r \frac{M}{m r^2} \frac{dr}{\sqrt{\frac{2}{M} ( E - U + \frac{M^2}{2 m r^2}) }}\end{aligned} \hspace{\stretch{1}}(2.1)

but this can be messy to solve, where we get elliptic integrals or worse, depending on the potential.

For the special case of the 3D problem where the potential has a 1/r form, this is what Prof. Poppitz called “super-integrable”. With 2N - 1 = 5 conserved quantities to be found, we’ve got one more. Here the form of that last conserved quantity is given, called the Runge-Lenz vector, and we verify that it is conserved.

Runge-Lenz vector

Given a potential

\begin{aligned}U = -\frac{\alpha}{r}\end{aligned} \hspace{\stretch{1}}(3.2)

and a Lagrangian

\begin{aligned}\mathcal{L} &= \frac{m \dot{r}^2}{2} + \frac{1}{{2}} \frac{M_z^2}{m r^2} - U \\ M_z &= m r^2 \dot{\phi}^2\end{aligned} \hspace{\stretch{1}}(3.3)

and writing the angular momentum as

\begin{aligned}\mathbf{M} = m \mathbf{r} \times \mathbf{v} \end{aligned} \hspace{\stretch{1}}(3.5)

the Runge-Lenz vector

\begin{aligned}\mathbf{A} = \mathbf{v} \times \mathbf{M} - \alpha \hat{\mathbf{r}},\end{aligned} \hspace{\stretch{1}}(3.6)

is a conserved quantity.

Verify the conservation assumption.

Let’s show that the conservation assumption is correct

\begin{aligned}\frac{d}{dt} \left( \mathbf{v} \times \mathbf{M} \right)=\frac{d{{ \mathbf{v}}}}{dt} \times \mathbf{M} + \mathbf{v} \times \not{{\frac{d{{\mathbf{M} }}}{dt}}}\end{aligned} \hspace{\stretch{1}}(3.7)

Here, we note that angular momentum conservation is really d\mathbf{M}/dt = 0, so we are left with only the acceleration term, which we can rewrite in terms of the Euler-Lagrange equation

\begin{aligned}\frac{d}{dt} \left( \mathbf{v} \times \mathbf{M} \right)&=-\frac{1}{{m}} \boldsymbol{\nabla} U \times M \\ &=-\frac{1}{{m}} \frac{\partial {U}}{\partial {r}} \hat{\mathbf{r}} \times M \\ &=-\frac{1}{{m}} \frac{\partial {U}}{\partial {r}} \hat{\mathbf{r}} \times (m \mathbf{r} \times \mathbf{v}) \\ &=- \frac{\partial {U}}{\partial {r}} \hat{\mathbf{r}} \times (\mathbf{r} \times \mathbf{v}) \end{aligned}

We can compute the double cross product

\begin{aligned}(\mathbf{a} \times (\mathbf{b} \times \mathbf{c}) )_i&=a_m b_r c_s \epsilon_{r s t} \epsilon_{m t i} \\ &=a_m b_r c_s \delta^{[rs]}_{i m} \\ &=a_m b_i c_m -a_m b_m c_i \end{aligned}

For

\begin{aligned}\mathbf{a} \times (\mathbf{b} \times \mathbf{c}) = (\mathbf{a} \cdot \mathbf{c}) \mathbf{b} -(\mathbf{a} \cdot \mathbf{b}) \mathbf{c}\end{aligned} \hspace{\stretch{1}}(3.8)

Plugging this we have

\begin{aligned}\frac{d}{dt} \left( \mathbf{v} \times \mathbf{M} \right) &= \frac{\partial U}{\partial r} \\ left( (\hat{\mathbf{r}} \cdot \mathbf{r}) \mathbf{v}-(\hat{\mathbf{r}} \cdot \mathbf{v}) \mathbf{r} \right) \\ &= \left( \frac{\alpha}{r^2} \right)\left( r \mathbf{v}-\frac{1}{r}(\mathbf{r} \cdot \mathbf{v}) \mathbf{r} \right) \\ &= \alpha\left( \frac{\mathbf{v}}{r}-\frac{(\mathbf{r} \cdot \mathbf{v}) \mathbf{r} }{r^3}\right) \\ \end{aligned}

Now let’s look at the other term. We’ll need the derivative of \hat{\mathbf{r}}

\begin{aligned}\frac{d{{\hat{\mathbf{r}}}}}{dt} &=\frac{d}{dt} \frac{\mathbf{r}}{r} \\ &=\frac{\mathbf{v}}{r} + \mathbf{r} \frac{d{{\frac{1}{{r}}}}}{dt} \\ &=\frac{\mathbf{v}}{r} - \frac{\mathbf{r}}{r^2} \frac{d{{r}}}{dt} \\ &=\frac{\mathbf{v}}{r} - \frac{\mathbf{r}}{r^2} \frac{d{{ \sqrt{\mathbf{r} \cdot \mathbf{r}}}}}{dt} \\ &=\frac{\mathbf{v}}{r} - \frac{\mathbf{r}}{r^2} \frac{\mathbf{v} \cdot \mathbf{r}}{\sqrt{\mathbf{r}^2}} \\ &=\frac{\mathbf{v}}{r} - \frac{\mathbf{r}}{r^3} \mathbf{v} \cdot \mathbf{r}\end{aligned}

Putting all the bits together we’ve now verified the conservation statement

\begin{aligned}\frac{d}{dt} \left(\mathbf{v} \times \mathbf{M} - \alpha \hat{\mathbf{r}}\right)=\alpha\left( \frac{\mathbf{v}}{r}-\frac{(\mathbf{r} \cdot \mathbf{v}) \mathbf{r} }{r^3}\right) -\alpha \left( \frac{\mathbf{v}}{r} - \frac{\mathbf{r}}{r^3} \mathbf{v} \cdot \mathbf{r} \right)= 0.\end{aligned} \hspace{\stretch{1}}(3.9)

With

\begin{aligned}\frac{d}{dt} \left( \mathbf{v} \times \mathbf{M} - \alpha \hat{\mathbf{r}} \right) = 0,\end{aligned} \hspace{\stretch{1}}(3.10)

our vector must be some constant vector. Let’s write this

\begin{aligned}\mathbf{v} \times \mathbf{M} - \alpha \hat{\mathbf{r}} = \alpha \mathbf{e},\end{aligned} \hspace{\stretch{1}}(3.11)

so that

\begin{aligned}\boxed{\mathbf{v} \times \mathbf{M} = \alpha \left(\mathbf{e} + \hat{\mathbf{r}} \right).}\end{aligned} \hspace{\stretch{1}}(3.12)

Dotting 3.12 with \mathbf{M} we find

\begin{aligned}\alpha \mathbf{M} \cdot \left(\mathbf{e} + \hat{\mathbf{r}} \right)&=\mathbf{M} \cdot (\mathbf{v} \times \mathbf{M}) \\ &= 0\end{aligned}

With \hat{\mathbf{r}} lying in the plane of the trajectory (perpendicular to \mathbf{M}), we must also have \mathbf{e} lying in the plane of the trajectory.

Now we can dot 3.12 with \mathbf{r} to find

\begin{aligned}\mathbf{r} \cdot (\mathbf{v} \times \mathbf{M}) &= \alpha \mathbf{r} \cdot \left(\mathbf{e} + \hat{\mathbf{r}} \right) \\  &= \alpha \cdot \left( r e \cos(\phi - \phi_0) + r \right) \\ \mathbf{M} \cdot (\mathbf{r} \times \mathbf{v}) &= \\ \mathbf{M} \cdot \frac{\mathbf{M}}{m} &= \\ \frac{\mathbf{M}^2}{m} &=\end{aligned}

This is

\begin{aligned}\frac{\mathbf{M}^2}{m} = \alpha r \left( 1 + e \cos(\phi - \phi_0) \right).\end{aligned} \hspace{\stretch{1}}(3.13)

This is a kind of curious implicit relationship, since \phi is also a function of r. Recall that the kinetic portion of our Lagrangian was

\begin{aligned}\frac{1}{{2}} m (\dot{r}^2 + r^2 \dot{\phi}^2 )\end{aligned} \hspace{\stretch{1}}(3.14)

so that our angular momentum was

\begin{aligned}M_\phi = \frac{\partial }{\partial {\dot{\phi}}} \left( \frac{1}{{2}} m r^2 \dot{\phi}^2 \right) = m r^2 \dot{\phi},\end{aligned} \hspace{\stretch{1}}(3.15)

with no \phi dependence in the Lagrangian we have

\begin{aligned}\frac{d}{dt} (m r^2 \dot{\phi}) = 0,\end{aligned} \hspace{\stretch{1}}(3.16)

or

\begin{aligned}\mathbf{M} = m r^2 \dot{\phi} \hat{\mathbf{z}} = \text{constant}\end{aligned} \hspace{\stretch{1}}(3.17)

Our dynamics are now fully specified, even if this not completely explicit

\begin{aligned}\boxed{\begin{aligned}r &= \frac{M^2}{m \alpha} \frac{1}{{1 + e \cos(\phi - \phi_0)}} \\ \frac{d\phi}{dt} &= \frac{M}{ m r^2}.\end{aligned}}\end{aligned} \hspace{\stretch{1}}(3.18)

What we can do is rearrange and separate variables

\begin{aligned}\frac{1}{{r^2}} = \frac{m^2 \alpha^2}{M^4} (1 + e \cos(\phi - \phi_0))^2 = \frac{m}{M} \frac{d\phi}{dt},\end{aligned} \hspace{\stretch{1}}(3.19)

to find

\begin{aligned}t - t_0 = \frac{M^3}{m \alpha^3} \int_{\phi_0}^\phi d\phi \frac{1}{{(1 + e \cos(\phi - \phi_0))^2}}=\frac{M^3}{m \alpha^3} \int_0^{\phi - \phi_0} du\frac{1}{{(1 + e \cos u)^2}}\end{aligned} \hspace{\stretch{1}}(3.20)

Now, at least \phi = \phi(t) is specified implicitly.

We can also use the first of these to determine the magnitude of the radial velocity

\begin{aligned}\frac{dr}{dt} &=-\frac{M^2}{m \alpha} \frac{1}{{(1 + e \cos(\phi - \phi_0))^2}} (-e \sin(\phi - \phi_0)) \frac{d\phi}{dt} \\ &=\frac{e M^2}{m \alpha} \frac{1}{{(1 + e \cos(\phi - \phi_0))^2}} \sin(\phi - \phi_0) \frac{M}{m r^2} \\ &=\frac{e M^3}{m^2 \alpha r^2} \frac{1}{{(1 + e \cos(\phi - \phi_0))^2}} \sin(\phi - \phi_0) \\ &=\frac{e M^3}{m^2 \alpha r^2} \left( \frac{ m r \alpha }{M^2} \right)^2 \sin(\phi - \phi_0) \\ &=\frac{e }{M } \sin(\phi - \phi_0),\end{aligned}

with this, we can also find the energy

\begin{aligned}E &= \dot{r}( m \dot{r}) + \dot{\phi} ( m r^2 \dot{\phi}) - \left( \frac{1}{{2}} m \dot{r}^2 + \frac{1}{{2}} m r^2 \dot{\phi}^2 - U \right) \\ &= \frac{1}{{2}} m \dot{r}^2 + \frac{1}{{2}} m r^2 \dot{\phi}^2 + U  \\ &= \frac{1}{{2}} m \dot{r}^2 + \frac{1}{{2}} m r^2 \dot{\phi}^2 - \frac{\alpha}{r} \\ &= \frac{1}{{2}} m \frac{e^2}{M^2} \sin^2(\phi - \phi_0) + \frac{1}{{2 m r^2 }} M^2 - \frac{\alpha}{r}.\end{aligned}

Or

\begin{aligned}E= \frac{m}{2 M^2} (\mathbf{e} \times \hat{\mathbf{r}})^2 + \frac{1}{{2 m r^2 }} M^2 - \frac{\alpha}{r}.\end{aligned} \hspace{\stretch{1}}(3.21)

Is this what was used in class to state the relation

\begin{aligned}e = \sqrt{1 + \frac{2 E M^2}{m \alpha^2}}.\end{aligned} \hspace{\stretch{1}}(3.22)

It’s not obvious exactly how that is obtained, but we can go back to 3.18 to eliminate the e^2 \sin^2 \Delta \phi term

\begin{aligned}E = \frac{1}{{2}} m \frac{1}{M^2} \left( e^2 - \left( \frac{M^2}{r m \alpha} - 1\right)^2 \right) + \frac{1}{{2 m r^2 }} M^2 - \frac{\alpha}{r}.\end{aligned} \hspace{\stretch{1}}(3.23)

Presumably this simplifies to the desired result (or there’s other errors made in that prevent that).

Posted in Math and Physics Learning. | Tagged: , , , , , , , | Leave a Comment »