Peeter Joot's Blog.

Math, physics, perl, and programming obscurity.

Posts Tagged ‘Mathematica’

« Previous Entries

Steady state velocity profile of stirred cup of non-bottomless coffee (continued).

Posted by peeterjoot on April 30, 2012

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

Solving for the flow below the stirring point.

Previously, we found the functional form for an azimuthal flow that has z-axis dependence. Attempting apply boundary conditions that including a cylindrical mixing device got us into trouble, since our no-slip condition would simultaneously require zero and non-zero velocity at the point of contact of the mixing interface and the bottom of the vessel (i.e. where the stir stick contacts the bottom of the cup). We can avoid this issue by constraining the mixing to only occur above the bottom of the cup, and then look at the flow that this induces below the mixing point.

Let’s attempt to solve the boundary value problem. We’ll position the mixing cylinder at a height d above the bottom of the cup, and set the radius of that inner cylinder to s as before, with the mixing occurring at an angular velocity of \Omega. We now want apply these boundary value constrains to the velocity function we’ve found for the steady state problem. Provided z < d, this will have the form

\begin{aligned}u(r, z, 0) = s \Omega \sum C_i J_1(\lambda_i r/R) \sinh(\lambda_i z/R)\end{aligned} \hspace{\stretch{1}}(1.1)

where \lambda_i are the zeros of the order one Bessel function J_1. Observe that our boundary value conditions of u(R, z, 0) = 0 and u(r, 0, 0) = 0 are automatically satisfied. Note that because of the u(R, z, 0) = 0 equality here, this form of solution is no good if we are mixing at the edge of the cup, so we require s \ne R.

Our remaining boundary value condition u(s, d, 0) = s \Omega, means that we have to solve

\begin{aligned}1 = \sum C_i J_1(\lambda_i s/R) \sinh(\lambda_i d/R)\end{aligned} \hspace{\stretch{1}}(1.2)

Having had the same sort of problem in our steady state bottomless coffee problem, we know how to solve this

\begin{aligned}C_i \sinh(\lambda_i d/R) = \frac{\int_0^1 w J_1 (\lambda_i w) dw}{\int_0^1 w J_1^2 (\lambda_i w) dw}.\end{aligned} \hspace{\stretch{1}}(1.3)

So we find that below the point of stirring (z = d), our steady state solution for the velocity should be

\begin{aligned}\mathbf{u}(r, z, 0) = s \Omega \hat{\boldsymbol{\phi}} \sum_{i=1}^\infty\frac{\int_0^1 w J_1 (\lambda_i w) dw}{\int_0^1 w J_1^2 (\lambda_i w) dw}J_1(\lambda_i r/R) \frac{\sinh(\lambda_i z/R)}{\sinh(\lambda_i d/R) }\end{aligned} \hspace{\stretch{1}}(1.4)

For a cup size of R = 5 \text{cm}, stir radius of s = 3 \text{cm}, stir depth d = 2 \text{cm}, and angular velocity \Omega = 2 \pi \text{rad}/\text{s}, we find

\begin{aligned}\begin{aligned}u(r, z, 0) &=27.1119 J_{1}(7.66341 r) \sinh (7.66341 z)-3.42819 J_{1}(14.0312 r) \sinh (14.0312 z) \\ &+4.97807 J_{1}(20.3469 r) \sinh (20.3469 z)-1.53542 J_{1}(26.6474 r) \sinh (26.6474 z) + \cdots\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.5)

We can plot this velocity field as in figure (1), but it’s hard to see the radial dependence.

Figure 1: Vector field plot of the velocity field below the stir depth.

The radial dependence of the magnitude of the velocity can be seen in figure (2), which plots u(r, d, 0). We see the zero velocity at the edges of the cup as expected, and once we hit the stir height, matching of stir velocity. An animation showing the variation of the radial velocity profile at various depths up to the stir height is available at

, and also below.

Figure 2: Radial velocity dependence at the stir height.

Observing this animation we see that the velocity is dominated by the first term in the Bessel series, essentially just scaled by the hyperbolic sine that multiplies it.

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

Flow between infinite moving inner cylinder and outer cylinder (3D visualizations).

Posted by peeterjoot on April 16, 2012

[Click here for a PDF of this post with nicer formatting.]

We previously looked at the problem of an infinite cylinder of radius R_1 is moving with velocity v parallel to its axis. It is places inside another cylinder of radius R_2. The axes of the two cylinders coincide. The fluid is incompressible, with viscosity \mu and density \rho, the flow is assumed to be stationary, and no external pressure gradient is applied, and found the solution to have the form

\begin{aligned}\boxed{w = \frac{G}{4 \mu} (R_2^2 - r^2) + \left(v - \frac{G}{4 \mu} (R_2^2 - R_1^2)\right)\frac{\ln r/R_2}{\ln R_1/R_2}.}\end{aligned} \hspace{\stretch{1}}(3.27)

This system with such simple geometry lends itself well to 3D visualization of the velocity profiles

An animation of this can be found at

. The Mathematica notebook that generated this is available at phy454/twoCylinders3D.cdf. The 3D visualization portion of that notebook has also been made available as part of the wolfram demonstrations project.

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

Couette flow of a viscous incompressible fluid (visualizing it).

Posted by peeterjoot on April 12, 2012

[Click here for a PDF of this post(and previous).]

In Couette flow of a viscous incompressible fluid, I looked at solving this problem. As a followup, here are some visualizations of the velocity field as the outer cylinder angular velocity is changed.

The Mathematica workbook that was used to generate the figures is also available for more direct play, and includes slider controls for the cylinder radii and angular velocities.

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

Compilation of class notes for phy454h1s, continuum mechanics (including all lectures)

Posted by peeterjoot on April 6, 2012

As mentioned before the midterm, I’ve got a compilation of class notes here:


https://sites.google.com/site/peeterjoot2/math2012/phy454.pdf

This has been updated now with the full set of lecture notes from class. Included in this document are:

  • my lecture notes, often auxmented by whatever I felt required to make them understandable to myself
  • some solved problems
  • some stuff I was playing with on the side
  • links to the various mathematica notebooks (some of which have some cool stuff, and can be viewed after download with the free wolfram CDF player).

All of the following, previously posted individually are included, and will no longer be maintained independently

Apr 4, 2012 Thermal stability.

Mar 31, 2012 Channel flow with step pressure gradient

Mar 29, 2012 Stability.

Mar 28, 2012 Asymptotic solutions of ill conditioned equations.

Mar 23, 2012 Boundary layers and stream functions. The Blasius problem.

Mar 23, 2012 Problem Set 3. Velocity scaling, non-dimensionalisation, boundary layers.

Mar 21, 2012 Scaling to setup to solve the boundary layer problem.

Mar 16, 2012 Impulsive flow. Boundary layers. Oscillatory driven flow.

Mar 14, 2012 Midterm reflection.

Mar 14, 2012 Boundary Layers. Time dependent flow.

Mar 9, 2012 More on surface tension and Reynold’s number.

Mar 7, 2012 Non-dimensionality and scaling.

Mar 6, 2012 Inclined flow without constant height assumption. Setting up the problem, but not actually solving it.

Mar 2, 2012 Problem Set 2. Rectilinear flow with pressure gradient. Two layer flow driven by moving surface.

Mar 2, 2012 Hydrostatics. Surface normals and tangent vectors.

Mar 1, 2012 Steady state inclined flow down a plane of two layers of incompressible viscous fluids.

Mar 1, 2012 Steady state inclined flow down a plane of two layers of incompressible viscous equal density fluids.

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: , , , | Leave a Comment »

Channel flow with step pressure gradient

Posted by peeterjoot on April 1, 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 and statement.

This is problem 2.5 from [1].

Viscous fluid is at rest in a two-dimensional channel between stationary rigid walls with y = \pm h. For t \ge 0 a constant pressure gradient P = -dp/dx is imposed. Show that u(y, t) satifies

\begin{aligned}\frac{\partial {u}}{\partial {t}} = \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2} + \frac{P}{\rho},\end{aligned} \hspace{\stretch{1}}(1.1)

and give suitable initial and boundary conditions. Find u(y, t) in form of a Fourier series, and show that the flow approximates to steady channel flow when t \gg h^2/\nu.

Solution

With only horizontal components to the flow, the Navier-Stokes equations for incompressible flow are

\begin{subequations}

\begin{aligned}\frac{\partial {u}}{\partial {t}} + u \frac{\partial {u}}{\partial {x}} = -\frac{1}{{\rho}} \frac{\partial {p}}{\partial {x}} + \nu \left( \frac{\partial^2 {{}}}{\partial {{x}}^2}+\frac{\partial^2 {{}}}{\partial {{y}}^2}\right)u\end{aligned} \hspace{\stretch{1}}(2.2a)

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

\end{subequations}

Substitution of 2.2b into 2.2a gives us

\begin{aligned}\frac{\partial {u}}{\partial {t}} = -\frac{1}{{\rho}} \frac{\partial {p}}{\partial {x}} + \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2}.\end{aligned} \hspace{\stretch{1}}(2.3)

Our equation to solve is therefore

\begin{aligned}\boxed{\frac{\partial {u}}{\partial {t}} = \Theta(t) \frac{P}{\rho} + \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2}.}\end{aligned} \hspace{\stretch{1}}(2.4)

This equation, for t < 0, allows for solutions

\begin{aligned}u = A y + B\end{aligned} \hspace{\stretch{1}}(2.5)

but the problem states that the fluid is at rest initially, so we don’t really have to solve anything (i.e. A = B = 0).

The no-slip conditions introduce boundary value conditions u(\pm h, t) = 0.

For t \ge 0 we have

\begin{aligned}\frac{\partial {u}}{\partial {t}} = \frac{P}{\rho} + \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2}.\end{aligned} \hspace{\stretch{1}}(2.6)

If we attempt separation of variables with u(y, t) = Y(y) T(t), our equation takes the form

\begin{aligned}T' Y = \frac{P}{\rho} + \nu T Y''.\end{aligned} \hspace{\stretch{1}}(2.7)

We see that the non-homogenous term prevents successful application of separation of variables. Let’s modify our problem by attempting to recast our equation into a homogenous form by adding a particular solution for the steady state flow problem. That problem was the solution of

\begin{aligned}\frac{\partial^2 {{}}}{\partial {{y}}^2} u_s(y, 0) = -\frac{P}{\rho \nu}\end{aligned} \hspace{\stretch{1}}(2.8)

which has solution

\begin{aligned}u_s(y, 0) = \frac{P}{2 \rho \nu} \left( h^2 - y^2 \right) + A y + B.\end{aligned} \hspace{\stretch{1}}(2.9)

The freedom to incoporate an h^2 constant into the equation as an integration constant has been employed, knowing that it will kill the y^2 contributions at y = \pm h to make the boundary condition matching easier. Our no-slip conditions give us

\begin{aligned}0 &= A h + B \\ 0 &= -A h + B.\end{aligned} \hspace{\stretch{1}}(2.10)

Adding this we have 2 B = 0, and subtracting gives us 2 A h = 0, so a specific solution that matches our required boundary value (and initial value) conditions is just the steady state channel flow solution we are familiar with

\begin{aligned}u_s(y, 0) = \frac{P}{2 \rho \nu} \left( h^2 - y^2 \right).\end{aligned} \hspace{\stretch{1}}(2.12)

Let’s now assume that our general solution has the form

\begin{aligned}u(y, t) = u_H(y, t) + u_s(y, 0).\end{aligned} \hspace{\stretch{1}}(2.13)

Applying the Navier-Stokes equation to this gives us

\begin{aligned}\frac{\partial {u_H}}{\partial {t}} = \frac{P}{\rho} + \nu \frac{\partial^2 {{u_H}}}{\partial {{y}}^2} + \nu \frac{\partial^2 {{u_s}}}{\partial {{y}}^2}.\end{aligned} \hspace{\stretch{1}}(2.14)

But from 2.8, we see that all we have left is a homogenous problem in u_H

\begin{aligned}\frac{\partial {u_H}}{\partial {t}} = \nu \frac{\partial^2 {{u_H}}}{\partial {{y}}^2},\end{aligned} \hspace{\stretch{1}}(2.15)

where our boundary value conditions are now given by

\begin{aligned}0 &= u_H(\pm h, t) + u_s(\pm h) \\ &= u_H(\pm h, t),\end{aligned}

and

\begin{aligned}0 &= u(y, 0) \\ &= u_H(y, 0) + \frac{P}{2 \rho \nu} \left( h^2 - y^2 \right),\end{aligned}

or
\begin{subequations}

\begin{aligned}u_H(\pm h, t) = 0\end{aligned} \hspace{\stretch{1}}(2.16a)

\begin{aligned}u_H(y, 0) = -\frac{P}{2 \rho \nu} \left( h^2 - y^2 \right).\end{aligned} \hspace{\stretch{1}}(2.16b)

\end{subequations}

Now we can apply separation of variables with u_H = T(t) Y(y), yielding

\begin{aligned}T' Y = \nu T Y'',\end{aligned} \hspace{\stretch{1}}(2.17)

or

\begin{aligned}\frac{T'}{T} = \nu \frac{Y''}{Y} = \text{constant} = - \nu \alpha^2.\end{aligned} \hspace{\stretch{1}}(2.18)

Here a positive constant \nu \alpha^2 has been used assuming that we want a solution that is damped with time.

Our solutions are

\begin{aligned}T &\propto e^{- \nu \alpha^2 t} \\ Y &= A \sin \alpha y + B \cos\alpha y,\end{aligned} \hspace{\stretch{1}}(2.19)

or

\begin{aligned}u_H(y, t) = \sum_\alpha e^{-\alpha^2 \nu t} \left( A_\alpha \sin \alpha y + B_\alpha \cos\alpha y \right).\end{aligned} \hspace{\stretch{1}}(2.21)

We have constraints on \alpha due to our boundary value conditions. For our sin terms to be solutions we require

\begin{aligned}\sin (\alpha (\pm h)) = \sin n \pi\end{aligned} \hspace{\stretch{1}}(2.22)

and for our cosine terms to be solutions we require

\begin{aligned}\cos (\alpha (\pm h)) = \cos \left( \frac{\pi}{2} + n \pi \right)\end{aligned} \hspace{\stretch{1}}(2.23)

\begin{aligned}\alpha &= \frac{2 n \pi}{2 h} \\ \alpha &= \frac{2 n + 1 \pi}{2 h}\end{aligned} \hspace{\stretch{1}}(2.24)

respectively.

Our homogeneous solution therefore takes the form

\begin{aligned}u_H(y, t) = C_0 + \sum_{m > 0} C_m e^{ -(m \pi/2h)^2 \nu t } \left\{\begin{array}{l l}\sin \left( \frac{ m \pi y }{2 h} \right) & \quad \mbox{latex m$ even} \\ \cos \left( \frac{ m \pi y }{2 h} \right) & \quad \mbox{m odd} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.26)$

Our undetermined constants should be provided by the boundary value constraint at t = 0 2.16b, leaving us to solve the Fourier problem

\begin{aligned}-\frac{P}{2 \mu} \left( h^2 - y^2 \right)=\sum_{m \ge 0} C_m \left\{\begin{array}{l l}\sin \left( \frac{ m \pi y }{2 h} \right) & \quad \mbox{latex m$ even} \\ \cos \left( \frac{ m \pi y }{2 h} \right) & \quad \mbox{m odd} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.27)$

Multiplying by a sine and integrating will clearly give zero (even times odd function over a symmetric interval). Let’s see if there’s any scaling required to select out the C_m term

\begin{aligned}\int_{-h}^h \cos \left( \frac{ m \pi y }{2 h} \right) \cos \left( \frac{ n \pi y }{2 h} \right) dy&=\frac{2h}{\pi} \int_{-h}^h \cos \left( \frac{ m \pi y }{2 h} \right) \cos \left( \frac{ n \pi y }{2 h} \right) \pi dy/2h \\ &=\frac{2h}{\pi} \int_{-\pi/2}^{\pi/2}\cos m x \cos n x dx \\ &=\frac{h}{\pi} \int_{-\pi/2}^{\pi/2}\left( \cos( (m - n) \pi/2 ) +\cos( (m + n) \pi/2 ) \right) dx\end{aligned}

Note that since m and n must be odd, m \pm n = 2 c for some integer c, so this integral is zero unless m = n (consider m = 2 a + 1, n = 2 b + 1). For the m = n term we have

\begin{aligned}\int_{-h}^h \cos \left( \frac{ m \pi y }{2 h} \right) \cos \left( \frac{ n \pi y }{2 h} \right) dy&=\frac{h}{\pi} \int_{-\pi/2}^{\pi/2}\left( 1+\cos( m \pi ) \right) dx \\ &=h\end{aligned}

Therefore, our constants C_m (for odd m) are given by

\begin{aligned}C_m &= -\frac{P h }{2 \mu} \int_{-h}^h\left( 1 - \left( \frac{y}{h}\right)^2 \right)\cos \left( \frac{ m \pi y }{2 h} \right) dy \\ &=-\frac{P h^2 }{2 \mu} \int_{-1}^1\left( 1 - x^2 \right)\cos \left( \frac{ m \pi x }{2} \right) x \\ \end{aligned}

With m = 2 n + 1, we have

\begin{aligned}C_{2 n + 1} = -\frac{16 P h^2 (-1)^n}{\mu \pi^3 (2 n + 1)^3}.\end{aligned} \hspace{\stretch{1}}(2.28)

For that calculation see: phy454/channelFlowWithStepPressureGradient.cdf.

Our complete solution is

\begin{aligned}u(y, t) =\frac{P h^2}{2 \mu} \left( 1 - \left( \frac{y}{h} \right)^2 \right)- \frac{16 P h^2 }{\mu \pi^3 }\sum_{n = 0}^\infty \frac{(-1)^n}{(2 n + 1)^3}e^{ -((2 n + 1) \pi/2h)^2 \nu t } \cos \left( \frac{ (2 n + 1) \pi y }{2 h} \right) .\end{aligned} \hspace{\stretch{1}}(2.29)

The largest of the damped exponentials above is the n = 0 term which is

\begin{aligned}e^{ - \pi^2 \nu t /h^2 },\end{aligned} \hspace{\stretch{1}}(2.30)

so if \nu t >> h^2 these terms all die off, leaving us with just the steady state.

Rather remarkably, this Fourier series is actually a very good fit even after only a single term. Using the viscosity and density of water, h = 1 \text{cm}, and P = 3 \times \mu_{\text{water}} \times (2 \text{cm}/{s})/ h^2 (parameterizing the pressure gradient by the average velocity it will induce), a plot of the parabola that we are fitting to and the difference of that from the first Fourier term is shown in figure (1).

Figure 1: Parabolic channel flow steady state, and difference from first Fourier term.

The higher order corrections are even smaller. Even the first order deviations from the parabola that we are fitting to is a correction on the scale of 1/100 of the height of the parabola. This is illustrated in figure (2) where the magnitude of the first 5 deviations from the steady state are plotted.

Figure 2: Difference from the steady state for the first five Fourier terms.

An animation of the time evolution above can be found in figure (3). If this animation is unavailable, it can also be found at

.

Figure 3: Time evolution of channel flow velocity profile after turning on a constant pressure gradient.

It’s also interesting to look at the very earliest part of the time evolution. Observe in figure (4)  (or

) the oscillatory phenomina. Could some of that be due to not running with enough Fourier terms in this early part of the evolution when more terms are probably significant?

Figure 4: Early time evolution of channel flow velocity profile after turning on a constant pressure gradient.

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. Problem Set 2. Two layer flow driven by moving surface. (revisited)

Posted by peeterjoot on March 26, 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.)]

Problem Q3 (revisited).

I’d produced the following sketches. For a higher viscosity bottom layer \mu^{(1)} > \mu^{(2)}, this should look something like figure (\ref{fig:continuumProblemSet2:continuumProblemSet2Fig4}) whereas for the higher viscosity on the top, these would be roughly flipped as in figure (\ref{fig:continuumProblemSet2:continuumProblemSet2Fig5}).

This superposition can be justified since we have no (\mathbf{u} \cdot \boldsymbol{\nabla})\mathbf{u} term in the Navier-Stokes equations for these systems.

Exact solutions.

The figures above are kind of rough. It’s not actually hard to solve the system above. After some simplification, I find with Mathematica the following solution

\begin{aligned}u^{(1)}(y) &= -\frac{\mu^{(2)} U -G h^2}{\mu^{(1)}+\mu^{(2)}}-\frac{y \left(G h^2 (\mu^{(2)} -\mu^{(1)}) +2 \mu^{(1)} \mu^{(2)} U \right)}{2 h \mu^{(1)} \left(\mu^{(1)}+\mu^{(2)}\right)}-\frac{G y^2}{2 \mu^{(1)}} \\ u^{(2)}(y) &= -\frac{\mu^{(2)} U -G h^2}{\mu^{(1)}+\mu^{(2)}}-\frac{y \left(G h^2 (\mu^{(2)} -\mu^{(1)}) +2 \mu^{(1)} \mu^{(2)} U \right)}{2 h \mu^{(2)} \left(\mu^{(1)}+\mu^{(2)}\right)}-\frac{G y^2}{2 \mu^{(2)}}.\end{aligned} \hspace{\stretch{1}}(3.49)

Should we wish a more exact plot for any specific values of the viscosities, we could plot exactly with software the vector field described by these velocities.

I suppose it is cheating to use Mathematica and then say that the solution is easy? To make amends for being lazy with my algebra, let’s show that it is easy to do manually too. I’ll do the same problem manually, but generalize it slightly. We can do this easily if we just be a bit marter with our integration constants. Let’s solve the problem for the upper and lower walls moving with velocity V_2 and V_1 respectively, and let the heights from the interface be h_2 and h_1 respectively.

We have the same set of differential equations to solve, but now let’s write our solution with the undetermined coefficients expressed as

\begin{aligned}u^{(2)} &= -\frac{G}{2 \mu^{(2)}}\left( h_2^2 - y^2 \right) + \frac{A_2}{\mu^{(2)}} (y - h_2) + B_2 \\ u^{(1)} &= -\frac{G}{2 \mu^{(1)}}\left( h_1^2 - y^2 \right) + \frac{A_1}{\mu^{(1)}} (y + h_1) + B_1.\end{aligned} \hspace{\stretch{1}}(3.55)

Now it’s super easy to match the boundary conditions at y = -h_1 and y = h_2(the lower and upper walls respectively). Clearly the integration constants B_1,B_2 are just the velocities. Matching the tangential component of the traction vectors at y = 0 we have

\begin{aligned}A_2 = A_1\end{aligned} \hspace{\stretch{1}}(3.55)

and matching velocities at y = 0 gives us

\begin{aligned}-\frac{G}{2 \mu^{(2)}}h_2^2 - \frac{A_2}{\mu^{(2)}} h_2 + V_2 = -\frac{G}{2 \mu^{(1)}}h_1^2 + \frac{A_1}{\mu^{(1)}} h_1 + V_1.\end{aligned} \hspace{\stretch{1}}(3.55)

This gives us

\begin{aligned}u^{(2)} &= \frac{G h_2^2}{2 \mu^{(2)}}\left(\frac{y^2}{h_2^2} - 1 \right) + A \frac{ h_2}{\mu^{(2)}} \left( \frac{y}{h_2} - 1 \right) + V_2 \\ u^{(1)} &= \frac{G h_1^2}{2 \mu^{(1)}}\left(\frac{y^2}{h_1^2} - 1 \right) + A \frac{ h_1}{\mu^{(1)}} \left( \frac{y}{h_1} + 1 \right) + V_1 \\ A &=\frac{V_2 - V_1+ \frac{G h_1^2}{2 \mu^{(1)}}-\frac{G h_2^2}{2 \mu^{(2)}}}{\frac{h_1}{\mu^{(1)}}+\frac{h_2}{\mu^{(2)}}}.\end{aligned} \hspace{\stretch{1}}(3.55)

Plotting this with sliders or animation in Mathematica is a fun way to explore visualizing this. The results vary widely depending on the various parameters. Such a Mathematica notebook can be found in
https://raw.github.com/peeterjoot/physicsplay/master/notes/phy454/problemSetIIQ3PlotWithManipulate.cdf
. A pair of cached animations with variation of the pressure gradient for v_1 = 0, h_1 = h_2, showing the superposition of the shear and channel flow solutions can be found here

  1. With \mu^{(1)} > \mu^{(2)} fluids with densities of mercury and water in the lower (1) and upper (2) layers respectively.
  2. With \mu^{(2)} > \mu^{(1)}, this is plotted fluids with densities of water and mercury in the lower (1) and upper (2) layers respectively.

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

Steady state inclined flow down a plane of two layers of incompressible viscous fluids.

Posted by peeterjoot on March 4, 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.

Here’s a generalization of one of the problems from section 2 of [1], itself a slight variation on what we did in class.

In the previous calculation we did the calculation for two incompressible fluids of the same densities flowing down an inclined plane. Now, let’s generalize this slightly, allowing for different densities.

I’m curious how much the air in the neighborhood of some flowing water gets dragged by that flow. It never occurred to me that this would occur, and I’d like to plug in some numbers and see what the results are. This should be something that can be modeled with two layers like this, one of fluid, one of air of a specific thickness (allowing pressure to vary due to the velocity gradient), and one final layer of air at a fixed pressure (atmospheric). I wouldn’t expect that problem to be much harder than this one, although it may end up being worthwhile to let a computer algebra system do some of the grunt work to solve all the resulting equations.

In the end, when we get to putting in some numbers for this problem, we can probably also get an idea how deep the region where the air gets dragged by the fluid can get.

Steady state inclined flow down a plane of two layers of incompressible viscous unequal density fluids.

Setup of the equations of motion for this system.

Our problem is illustrated in figure (\ref{fig:twoLayerInclinedFlowDifferentDensities:twoLayerInclinedFlowDifferentDensitiesFig1}) with a plane set at angle \alpha, fluid depths of h^{(1)} and h^{(2)} respectively, and viscosities \mu^{(1)} and \mu^{(2)}. We’ll write H = h^{(1)} + h^{(2)}. We have a pair of Navier-Stokes equations to solve

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{twoLayerInclinedFlowFig1}
\caption{Two fluids layers in inclined flow.}
\end{figure}

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

Our steady state and incompressibility constraints break this into a few independent equations

\begin{aligned}\rho^{(i)} \frac{\partial {\mathbf{u}^{(i)}}}{\partial {t}} &= 0 \\ \boldsymbol{\nabla} \cdot \mathbf{u}^{(i)} &= 0 \\ \rho^{(i)} (\mathbf{u}^{(i)} \cdot \boldsymbol{\nabla}) \mathbf{u}^{(i)} &= - \boldsymbol{\nabla} p^{(i)} + \mu^{(i)} \boldsymbol{\nabla}^2 \mathbf{u}^{(i)} + \rho^{(i)} \mathbf{g}\end{aligned} \hspace{\stretch{1}}(2.2)

Let’s require no components of the flows in the y, or z directions initially. As the equivalent of Newton’s law for fluid flows, conservation of linear momentum requires that for our steady state problem we have u_y = u_z = 0 for the flows in both fluid layers.

FIXME: should go back and think about momentum conservation in the context of fluids. I’m now so used to thinking about this as a symmetry issue from a Lagrangian context. With no Lagrangian here, what would be a robust way treat this in fluids? What is the Lagrangian for a fluid mechanics system? I’d imagine that it would be possible to set up a field Lagrangian with velocity fields.

Our problem is now reduced to a problem in four quantities (two velocities and two pressures). With \mathbf{u}^{(i)} = (u^{(i)}, 0, 0) we can restate Navier-Stokes in coordinate form as

\begin{subequations}

\begin{aligned}\partial_x u^{(i)}(x,y,z) = 0\end{aligned} \hspace{\stretch{1}}(2.5a)

\begin{aligned}\rho^{(i)} (u^{(i)} \partial_x u^{(i)} = - \partial_x p^{(i)} + \mu^{(i)} (\partial_{xx} + \partial_{yy} + \partial_{zz}) u^{(i)} + \rho^{(i)} g \sin\alpha \end{aligned} \hspace{\stretch{1}}(2.5b)

\begin{aligned}0 = - \partial_y p^{(i)} - \rho^{(i)} g \cos\alpha \end{aligned} \hspace{\stretch{1}}(2.5c)

\begin{aligned}0 = - \partial_z p^{(i)} \end{aligned} \hspace{\stretch{1}}(2.5d)

\end{subequations}

In order to solve this, we have eight simultaneous non-linear PDEs, four unknown functions, plus boundary conditions!

What are the boundary conditions? One is the “no-slip” condition, the experimental observation that velocities match at the interfaces. So we should have zero velocity for the fluid lying against the plane, and velocity matching between the fluids. The air above the fluid will also be flowing along at the rate of the uppermost portion of the top layer, but we’ll neglect that effect (i.e. considering two layers of equal density and not three, with one having a separate density). We also have matching of the traction vectors at the interfaces.

Writing this, it occurred to me that I didn’t fully understand what motivated the traction vector matching boundary value condition. Talking to our Prof about this, the matching of the traction vectors at any point can be thought of as an observational issue, but this is also a force balance issue. There is an induced velocity in the direction of the traction vector at any given point. For example, when we have unidirectional flow, we must have no normal component of the traction vector, and only a tangential component, because we have only the tangential flow. It is probably reasonable to think about this roughly as the equivalent of matching both acceleration and velocity at the boundary, but because densities and viscosities vary, we have to match the traction vectors and not the acceleration itself.

Before continuing to solve our Navier-Stokes equations let’s express the condition that the tangential component of the traction vectors match algebraically.

Dropping indexes temporarily, for the normal to the surface \hat{\mathbf{n}} = (n_1, n_2, n_3) = (0, 1, 0) we want to compute

\begin{aligned}T_1 &= \sigma_{1 k} n_k \\ &= \sigma_{1 k} \delta_{2 k} \\ &= \sigma_{1 2} \\ &= -p \not{{\delta_{1 2}}} + 2 \mu e^{1 2} \\ &= \mu \left( \frac{\partial {u_1}}{\partial {y}} + \not{{\frac{\partial {u_2}}{\partial {x}}}} \right) \end{aligned}

So the tangential component of the traction vector is

\begin{aligned}\mathbf{T}^{(i)} = \mu^{(i)} \frac{\partial {u^{(i)}}}{\partial {y}} \hat{\mathbf{x}}.\end{aligned} \hspace{\stretch{1}}(2.6)

As noted above, this is in fact, the only component of the traction vector, since we don’t have any non-horizontal flow.

Our boundary value conditions, what we need in addition to the Navier-Stokes equations of 2.5, to solve our problem, are the matching at any interface of the following conditions

\begin{subequations}

\begin{aligned}u^{(i)} = u^{(j)} \end{aligned} \hspace{\stretch{1}}(2.7a)

\begin{aligned}p^{(i)} = p^{(j)} \end{aligned} \hspace{\stretch{1}}(2.7b)

\begin{aligned}\mu^{(i)} \frac{\partial {u^{(i)}}}{\partial {y}} = \mu^{(i)} \frac{\partial {u^{(j)}}}{\partial {y}}.\end{aligned} \hspace{\stretch{1}}(2.7c)

\end{subequations}

There are actually three interfaces to consider, that of the lower layer liquid with the inclined plane, the interface between the two fluid layers, and the interface between the upper layer fluid and the air above it.

Solving our equations of motion.

Starting with the simplest, the z-coordinate equation, of Navier-Stokes 2.5d, we can conclude that each of the pressures is not a function of z, so that we have

\begin{aligned}p^{(i)} = p^{(i)}(x, y).\end{aligned} \hspace{\stretch{1}}(2.8)

Using this, we can integrate our y-coordinate Navier-Stokes equation 2.5c, to find

\begin{aligned}p^{(i)} = - \rho^{(i)} g y \cos\alpha + f^{(i)}(x)\end{aligned} \hspace{\stretch{1}}(2.9)

At this point we can introduce the first boundary value constraint, that the pressures must match at the interfaces. In particular, on the upper surface, where we have atmospheric pressure p_A our pressure is

\begin{aligned}p^{(2)}(H) = - \rho^{(2)} g H \cos\alpha + f^{(2)}(x) = p_A,\end{aligned} \hspace{\stretch{1}}(2.10)

so f^{(2)} is constant with value

\begin{aligned}f^{(2)}(x) = p_A + \rho^{(2)} g H \cos\alpha,\end{aligned} \hspace{\stretch{1}}(2.11)

which fully determines the density of the upper surface

\begin{aligned}p^{(2)}(y) = \rho^{(2)} g \cos\alpha (H - y) + p_A.\end{aligned} \hspace{\stretch{1}}(2.12)

Matching the pressure between the two layers of fluids we have

\begin{aligned}p^{(1)}(h_1) &= - \rho^{(1)} g h_1 \cos\alpha + f^{(1)}(x) \\ &= p^{(2)}(h_1) \\ &= \rho^{(2)} g \cos\alpha (H - h_1) + p_A \\ &= \rho^{(2)} g h_2 \cos\alpha + p_A,\end{aligned}

so that our undetermined function f^{(1)}(x)

\begin{aligned}f^{(1)}(x) = \left(\rho^{(1)} h_1 + \rho^{(2)} h_2 \right) g \cos\alpha + p_A.\end{aligned} \hspace{\stretch{1}}(2.13)

With the densities not equal, we no longer find that the pressure is dependent only on the total height y, independent of the velocities and viscosities

\begin{aligned}p^{(1)}(y) = g \cos\alpha \left(\rho^{(1)} (h_1 - y) + \rho^{(2)} h_2\right)+ p_A.\end{aligned} \hspace{\stretch{1}}(2.14)

However, this is still a fairly satisfying result. The pressure on the bottom layer is the total pressure due to the layer above it (the contribution due to the total height h_2 of that layer of the fluid). To that we add the pressure at our specific height, a linear function of the difference from the interface above it. Specified piecewise our pressure is now fully determined

\begin{aligned}\boxed{p(y) =\left\{\begin{array}{l l}g \cos\alpha \left(\rho^{(1)} (h_1 - y) + \rho^{(2)} h_2\right)+ p_A & \quad \mbox{latex y H$}\end{array}\right.}\end{aligned} \hspace{\stretch{1}}(2.15)$

Observe that we have the usual \rho g h form in all the terms of the pressure above, just scaled by the cosine of the angle since only a portion of the gravitational force is pushing normally on the fluids.

Having solved for the pressure, we are now set to return to the remaining Navier-Stokes equations 2.5a, and 2.5b for this system. From 2.5a we see that the non-linear term on the LHS of 2.5b is killed and also see that our velocities can only be functions of y and z

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

While more general solutions can likely be found, we will limit ourselves to looking only for solutions that are functions of y. From our solution to the pressure part of the problem p^{(i)} = p^{(i)}(y), we also see that the pressure term \partial_x p^{(i)} of 2.5b is killed. We are left with just

\begin{aligned}0 &= \mu^{(i)} (\partial_{xx} + \partial_{yy} + \partial_{zz}) u^{(i)} + \rho^{(i)} g \sin\alpha \\ &= \mu^{(i)} \partial_{yy} u^{(i)} + \rho^{(i)} g \sin\alpha \\ &= \mu^{(i)} \frac{d^2}{dy^2} u^{(i)} + \rho^{(i)} g \sin\alpha \end{aligned}

This is directly integrable, and we find for the velocities and traction vectors respectively

\begin{subequations}

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

\begin{aligned}T_x^{(i)} = \mu^{(i)} \frac{\partial {u^{(i)}}}{\partial {y}} = - \rho^{(i)} g y \sin\alpha + \mu^{(i)} A^{(i)} \end{aligned} \hspace{\stretch{1}}(2.17b)

\end{subequations}

The boundary conditions left to exploit are

\begin{aligned}u^{(1)}(0) &= 0 \\ u^{(1)}(h_1) &= u^{(2)}(h_1) \\ T_x^{(1)}(h_1) &= T_x^{(2)}(h_1) \\ T_x^{(2)}(H) &= 0,\end{aligned} \hspace{\stretch{1}}(2.18)

The first is the no-slip condition with the plane. The last is an approximation that assumes the liquid isn’t producing a measurable force on the air above it. The other two are for the interfaces between the two fluids.

From u^{(1)}(0) = 0 we see immediately that we have B^{(1)} = 0. From the traction vector equality in the atmosphere, we have

\begin{aligned}0 = - \rho^{(2)} g H \sin\alpha + \mu^{(2)} A^{(2)},\end{aligned} \hspace{\stretch{1}}(2.22)

or

\begin{aligned}A^{(2)} = \frac{\rho^{(2)} g H \sin\alpha }{ \mu^{(2)} }.\end{aligned} \hspace{\stretch{1}}(2.23)

These reduce the problem to solving for two last integration constants, where our velocities are

\begin{aligned}u^{(1)} &= -\frac{\rho^{(1)} g \sin\alpha }{2 \mu^{(1)}} y^2 + A^{(1)} y \\ u^{(2)} &= \frac{\rho^{(2)} g \sin\alpha }{2 \mu^{(2)}} \left( 2 H y -y^2 \right) + B^{(2)}.\end{aligned} \hspace{\stretch{1}}(2.24)

and our traction vectors are

\begin{aligned}T_x^{(1)} &= - \rho^{(1)} g y \sin\alpha + \mu^{(1)} A^{(1)} \\ T_x^{(2)} &= \rho^{(2)} g \sin\alpha \left( H - y \right).\end{aligned} \hspace{\stretch{1}}(2.26)

Matching both at the interface (y = h_1) gives us

\begin{aligned}-\frac{\rho^{(1)} g \sin\alpha }{2 \mu^{(1)}} h_1^2 + A^{(1)} h_1 &= \frac{\rho^{(2)} g \sin\alpha }{2 \mu^{(2)}} h_1 \left( 2 h_2 + h_1 \right) + B^{(2)} \\ - \rho^{(1)} g h_1 \sin\alpha + \mu^{(1)} A^{(1)} &= \rho^{(2)} g h_2 \sin\alpha \end{aligned} \hspace{\stretch{1}}(2.28)

We find

\begin{aligned}A^{(1)} = \frac{1}{{ \mu^{(1)} }}(\rho^{(1)} h_1 + \rho^{(2)} h_2 )g \sin\alpha \end{aligned} \hspace{\stretch{1}}(2.30)

Let’s substitute this back for our first fluid velocity

\begin{aligned}u^{(1)} &= -\frac{\rho^{(1)} g \sin\alpha }{2 \mu^{(1)}} y^2 + \frac{y}{ \mu^{(1)} }(\rho^{(1)} h_1 + \rho^{(2)} h_2 )g \sin\alpha \\ &=g \sin\alpha \left( -\frac{\rho^{(1)} }{2 \mu^{(1)}} y^2 + \frac{y}{ \mu^{(1)} }(\rho^{(1)} h_1 + \rho^{(2)} h_2 )\right) \\ &=\frac{g y \sin\alpha}{2 \mu^{(1)}} \left( \rho^{(1) } (2 h_1 - y) +2 \rho^{(2)} h_2 \right) \\ \end{aligned}

As a check we see this is consistent with the previous calculation when \rho^{(1)} = \rho^{(2)}. For our final integration constant we now find

\begin{aligned}B^{(2)} &=\frac{g h_1 \sin\alpha}{2 \mu^{(1)}} \left( \rho^{(1) } h_1 +2 \rho^{(2)} h_2 \right) -\frac{\rho^{(2)} g \sin\alpha }{2 \mu^{(2)}} h_1 \left( 2 h_2 + h_1 \right) \\ &=\frac{g h_1 \sin\alpha}{2}\left(\frac{1}{{\mu^{(1)}}} \left( \rho^{(1) } h_1 +2 \rho^{(2)} h_2 \right) -\frac{\rho^{(2)} }{\mu^{(2)}} \left( 2 h_2 + h_1 \right) \right) \\ &=\frac{g h_1 \sin\alpha}{2 \mu^{(2)}}\left(\frac{\mu^{(2)}}{\mu^{(1)}} \left( \rho^{(1) } h_1 +2 \rho^{(2)} h_2 \right) - \rho^{(2)} \left( 2 h_2 + h_1 \right) \right) \\ \end{aligned}

So, finally, we have for the velocities

\begin{subequations}

\begin{aligned}u^{(1)}(y) = \frac{g y \sin\alpha}{2 \mu^{(1)}} \left( \rho^{(1) } (2 h_1 - y) +2 \rho^{(2)} h_2 \right) \end{aligned} \hspace{\stretch{1}}(2.31a)

\begin{aligned}u^{(2)}(y) = \frac{g \sin\alpha }{2 \mu^{(2)}} \left(\rho^{(2)} \left( 2 H y -y^2 \right) + h_1 \left(\frac{\mu^{(2)}}{\mu^{(1)}} \left( \rho^{(1) } h_1 +2 \rho^{(2)} h_2 \right) - \rho^{(2)} \left( 2 h_2 + h_1 \right) \right) \right)\end{aligned} \hspace{\stretch{1}}(2.31b)

\end{subequations}

The final result looks reasonable. If the viscosities and densities are equal then we have the same velocity profile in both layers. That makes sense given the equal densities, since there would really be nothing that would then distinguish the two layers.

Numerical application.

To try this out numerically, I’ve created a Mathematica worksheet. The results are fairly suprising, showing either an error in this calculation, an error programming the worksheet, or the folly of even considering a steady state flow of this form for anything that is not extremely viscous. Some validation is required to see what’s up.

References

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

Posted in Math and Physics Learning. | Tagged: , , , , , , | 4 Comments »

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 »

PHY456H1S Continuum mechanics. Problem Set 1. Stress, Strain, Traction vector. Force free equilibrium.

Posted by peeterjoot on February 9, 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.

This problem set is as yet ungraded.

Problem Q1.

Statement

For the stress tensor

\begin{aligned}\sigma =\begin{bmatrix}6 & 0 & 2 \\ 0 & 1 & 1 \\ 2 & 1 & 3\end{bmatrix}\text{M Pa}\end{aligned} \hspace{\stretch{1}}(2.1)

Find the corresponding strain tensor, assuming an isotropic solid with Young’s modulus E = 200 \times 10^9 \text{N}/\text{m}^2 and Poisson’s ration \nu = 0.35.

Solution

We need to express the relation between stress and strain in terms of Young’s modulus and Poisson’s ratio. In terms of Lam\’e parameters our model for the relations between stress and strain for an isotropic solid was given as

\begin{aligned}\sigma_{ij} = \lambda e_{kk} \delta_{ij} + 2 \mu e_{ij}.\end{aligned} \hspace{\stretch{1}}(2.2)

Computing the trace

\begin{aligned}\sigma_{kk} = (3 \lambda + 2 \mu) e_{kk},\end{aligned} \hspace{\stretch{1}}(2.3)

allows us to invert the relationship

\begin{aligned}2 \mu e_{ij} = \sigma_{ij} - \lambda \frac{\sigma_{kk}}{3 \lambda + 2 \mu} \delta_{ij}.\end{aligned} \hspace{\stretch{1}}(2.4)

In terms of Poisson’s ratio \nu and Young’s modulus E, our Lam\’e parameters were found to be

\begin{aligned}\lambda &= \frac{ E \nu }{(1 - 2 \nu)(1 + \nu)} \\ \mu &= \frac{E}{2(1 + \nu)},\end{aligned} \hspace{\stretch{1}}(2.5)

and

\begin{aligned}3 \lambda + 2 \mu&= \frac{ 3 E \nu }{(1 - 2 \nu)(1 + \nu)} + \frac{E}{1 + \nu} \\ &= \frac{E}{1 + \nu} \left( \frac{3 \nu}{1 - 2 \nu} + 1\right) \\ &= \frac{E}{1 + \nu} \frac{1 + \nu}{1 - 2 \nu} \\ &= \frac{E}{1 - 2 \nu}.\end{aligned}

Our stress strain model for the relationship for an isotropic solid becomes
we find

\begin{aligned}\frac{E}{1 + \nu} e_{ij}&=\sigma_{ij}-\frac{ E \nu }{(1 - 2 \nu)(1 + \nu)} \frac{1 - 2 \nu}{E}\sigma_{kk} \delta_{ij} \\ &=\sigma_{ij}-\frac{ \nu }{1 + \nu}\sigma_{kk} \delta_{ij} \\ \end{aligned}

or

\begin{aligned}e_{ij}=\frac{1}{{E}}\left((1 + \nu)\sigma_{ij}-\nu\sigma_{kk} \delta_{ij}\right).\end{aligned} \hspace{\stretch{1}}(2.7)

As a sanity check note that this matches (5.12) of [1], although they use a notation of \sigma instead of \nu for Poisson’s ratio. We are now ready to tackle the problem. First we need the trace of the stress tensor

\begin{aligned}\sigma_{kk} = (6 + 1 + 3) \text{M Pa} = 10 \text{M Pa},\end{aligned} \hspace{\stretch{1}}(2.8)

\begin{aligned}e_{ij}&=\frac{1}{{E}}\left((1 + \nu)\begin{bmatrix}6 & 0 & 2 \\ 0 & 1 & 1 \\ 2 & 1 & 3\end{bmatrix}-10 \nu\begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix}\right)\text{M Pa} \\ &=\frac{1}{{E}}\left(\begin{bmatrix}6 & 0 & 2 \\ 0 & 1 & 1 \\ 2 & 1 & 3\end{bmatrix}+ 0.35\begin{bmatrix}-4 & 0 & 2 \\ 0 & -9 & 1 \\ 2 & 1 & -7\end{bmatrix}\right)\text{M Pa} \\ &=\frac{1}{{2 \times 10^{5}}}\left(\begin{bmatrix}6 & 0 & 2 \\ 0 & 1 & 1 \\ 2 & 1 & 3\end{bmatrix}+ 0.35\begin{bmatrix}-4 & 0 & 2 \\ 0 & -9 & 1 \\ 2 & 1 & -7\end{bmatrix}\right)\end{aligned}

Expanding out the last bits of arithmetic the strain tensor is found to have the form

\begin{aligned}e_{ij}=\begin{bmatrix} 23 & 0 & 13.5 \\ 0 & -10.75 & 6.75 \\ 13.5 & 6.75 & 2.75\end{bmatrix} 10^{-6}.\end{aligned} \hspace{\stretch{1}}(2.9)

Note that this is dimensionless, unlike the stress.

Problem Q2.

Statement

Small displacement field in a material is given by

\begin{aligned}e_1 &= 2 x_1 x_2 \\ e_2 &= x_3^2 \\ e_3 &= x_1^2 - x_3\end{aligned} \hspace{\stretch{1}}(3.10)

Find

\begin{enumerate}
\item the infinitesimal strain tensor e_{ij},
\item the principal strains and the corresponding principal axes at (x_1, x_2, x_3) = (1, 2, 4),
\item Is the body under compression or expansion?
\end{enumerate}

Solution. infinitesimal strain tensor e_{ij}

Diving right in, we have

\begin{aligned}e_{11}&= \frac{\partial {e_1}}{\partial {x_1}} \\ &= \frac{\partial {}}{\partial {x_1}}2 x_1 x_2 \\ &= 2 x_2\end{aligned}

\begin{aligned}e_{22}&= \frac{\partial {e_2}}{\partial {x_2}} \\ &= \frac{\partial {}}{\partial {x_2}} x_3^2 \\ &= 0\end{aligned}

\begin{aligned}e_{33}&= \frac{\partial {e_3}}{\partial {x_3}} \\ &= \frac{\partial {}}{\partial {x_3}} ( x_1^2 - x_3 ) \\ &= -1\end{aligned}

\begin{aligned}e_{12}&=\frac{1}{{2}} \left(\frac{\partial {e_2}}{\partial {x_1}}+\frac{\partial {e_1}}{\partial {x_2}}\right) \\ &=\frac{1}{{2}}\left(\not{{\frac{\partial {}}{\partial {x_1}} x_3^2 }}+\frac{\partial {}}{\partial {x_2}} 2 x_1 x_2\right) \\ &=x_1\end{aligned}

\begin{aligned}e_{23}&=\frac{1}{{2}} \left(\frac{\partial {e_3}}{\partial {x_2}}+\frac{\partial {e_2}}{\partial {x_3}}\right) \\ &=\frac{1}{{2}}\left(\not{{\frac{\partial {}}{\partial {x_2}} (x_1^2 - x_3 )}}+\frac{\partial {}}{\partial {x_3}} x_3^2\right) \\ &=x_3\end{aligned}

\begin{aligned}e_{31}&=\frac{1}{{2}} \left(\frac{\partial {e_1}}{\partial {x_3}}+\frac{\partial {e_3}}{\partial {x_1}}\right) \\ &=\frac{1}{{2}}\left(\not{{\frac{\partial {}}{\partial {x_3}} 2 x_1 x_2 }}+\frac{\partial {}}{\partial {x_1}} (x_1^2 - x_3 )\right) \\ &=x_1\end{aligned}

In matrix form we have

\begin{aligned}\mathbf{e} =\begin{bmatrix}2 x_2 & x_1 & x_1 \\ x_1 & 0 & x_3 \\ x_1 & x_3 & -1 \\ \end{bmatrix}\end{aligned} \hspace{\stretch{1}}(3.13)

Solution. principle strains and axes

At the point (1, 2, 4) the strain tensor has the value

\begin{aligned}\mathbf{e} =\begin{bmatrix}4 & 1 & 1 \\ 1 & 0 & 4 \\ 1 & 4 & -1\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(3.14)

We wish to diagonalize this, solving the characteristic equation for the eigenvalues \lambda

\begin{aligned}0 &=\begin{vmatrix}4 -\lambda & 1 & 1 \\ 1 & -\lambda & 4 \\ 1 & 4 & -1 -\lambda\end{vmatrix} \\ &=(4 -\lambda )\begin{vmatrix} -\lambda & 4 \\ 4 & -1 -\lambda\end{vmatrix}-\begin{vmatrix}1 & 1 \\ 4 & -1 -\lambda\end{vmatrix}+\begin{vmatrix}1 & 1 \\ -\lambda & 4 \\ \end{vmatrix} \\ &=(4 - \lambda)(\lambda^2 + \lambda - 16)-(-1 -\lambda - 4)+(4 + \lambda) \\ \end{aligned}

We find the characteristic equation to be

\begin{aligned}0 = -\lambda^3 + 3 \lambda^2 + 22\lambda - 55.\end{aligned} \hspace{\stretch{1}}(3.15)

This doesn’t appear to lend itself easily to manual solution (there are no obvious roots to factor out). As expected, since the matrix is symmetric, a plot (\ref{fig:continuumL8:continuumProblemSet1Q2fig1}) shows that all our roots are real

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumProblemSet1Q2fig1}
\caption{Q2. Characteristic equation.}
\end{figure}

Numerically, we determine these roots to be

\begin{aligned}\{5.19684, -4.53206, 2.33522\}\end{aligned} \hspace{\stretch{1}}(3.16)

with the corresponding basis (orthonormal eigenvectors), the principle axes are

\begin{aligned}\left\{\hat{\mathbf{p}}_1,\hat{\mathbf{p}}_2,\hat{\mathbf{p}}_3\right\}=\left\{\begin{bmatrix}0.76291 \\ 0.480082 \\ 0.433001\end{bmatrix},\begin{bmatrix}-0.010606 \\ -0.660372 \\ 0.750863\end{bmatrix},\begin{bmatrix}-0.646418 \\ 0.577433 \\ 0.498713\end{bmatrix}\right\}.\end{aligned} \hspace{\stretch{1}}(3.17)

Solution. Is body under compression or expansion?

To consider this question, suppose that as in the previous part, we determine a basis for which our strain tensor e_{ij} = p_i \delta_{ij} is diagonal with respect to that basis at a given point \mathbf{x}_0. We can then simplify the form of the stress tensor at that point in the object

\begin{aligned}\sigma_{ij}&=\frac{E}{1 + \nu} \left(e_{ij} + \frac{\nu}{1 - 2 \nu} e_{mm} \delta_{ij}\right) \\ &=\frac{E}{1 + \nu} \left(p_i + \frac{\nu}{1 - 2 \nu} e_{mm}\right)\delta_{ij}.\end{aligned}

We see that the stress tensor at this point is also necessarily diagonal if the strain is diagonal in that basis (with the implicit assumption here that we are talking about an isotropic material). Noting that the Poisson ratio is bounded according to

\begin{aligned}-1 \le \nu \le \frac{1}{{2}},\end{aligned} \hspace{\stretch{1}}(3.18)

so if our trace is positive (as it is in this problem for all points x_2 > 1/2), then any positive principle strain value will result in a positive stress along that direction). For example at the point (1,2,4) of the previous part of this problem (for which x_2 > 1/2), we have

\begin{aligned}\sigma_{ij}=\frac{E}{1 + \nu}\begin{bmatrix}5.19684+ \frac{3 \nu}{1 - 2 \nu} & 0 & 0 \\ 0 & -4.53206+ \frac{3 \nu}{1 - 2 \nu} & 0 \\ 0 & 0 & 2.33522+ \frac{3 \nu}{1 - 2 \nu}\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(3.19)

We see that at this point the (1,1) and (3,3) components of stress is positive (expansion in those directions) regardless of the material, and provided that

\begin{aligned}\frac{3 \nu}{1 - 2 \nu} > 4.53206\end{aligned} \hspace{\stretch{1}}(3.20)

(i.e. \nu > 0.375664) the material is under expansion in all directions. For \nu < 0.375664 the material at that point is expanding in the \hat{\mathbf{p}}_1 and \hat{\mathbf{p}}_3 directions, but under compression in the \hat{\mathbf{p}}_2 directions.

(save to disk and run with either Mathematica or the free Wolfram CDF player ( http://www.wolfram.com/cdf-player/ ) )

For a Mathematica notebook that visualizes this part of this problem see
https://raw.github.com/peeterjoot/physicsplay/master/notes/phy454/continuumProblemSet1Q2animated.cdf
. This animates the stress tensor associated with the problem, for different points (x,y,z) and values of Poisson’s ratio \nu, with Mathematica manipulate sliders available to alter these (as well as a zoom control to scale the graphic, keeping the orientation and scale fixed with any variation of the other parameters). This generalizes the solution of the problem (assuming I got it right for the specific (1,2,4) point of the problem). The vectors are the orthonormal eigenvectors of the tensor, scaled by the magnitude of the eigenvectors of the stress tensor (also diagonal in the basis of the diagonalized strain tensor at the point in question). For those directions that are under expansive stress, I’ve colored the vectors blue, and for compressive directions, I’ve colored the vectors red.

This requires either a Mathematica client or the free Wolfram CDF player, either of which can run the notebook after it is saved to your computer’s hard drive.

Problem Q3.

Statement

The stress tensor at a point has components given by

\begin{aligned}\sigma =\begin{bmatrix}1 & -2 & 2 \\ -2 & 3 & 1 \\ 2 & 1 & -1\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(4.21)

Find the traction vector across an area normal to the unit vector

\begin{aligned}\hat{\mathbf{n}} = ( \sqrt{2} \mathbf{e}_1 - \mathbf{e}_2 + \mathbf{e}_3)/2\end{aligned} \hspace{\stretch{1}}(4.22)

Can you construct a tangent vector \boldsymbol{\tau} on this plane by inspection? What are the components of the force per unit area along the normal \hat{\mathbf{n}} and tangent \boldsymbol{\tau} on that surface? (hint: projection of the traction vector.)

Solution

The traction vector, the force per unit volume that holds a body in equilibrium, in coordinate form was

\begin{aligned}P_i = \sigma_{ik} n_k\end{aligned} \hspace{\stretch{1}}(4.23)

where n_k was the coordinates of the normal to the surface with area df_k. In matrix form, this is just

\begin{aligned}\mathbf{P} = \sigma \hat{\mathbf{n}},\end{aligned} \hspace{\stretch{1}}(4.24)

so our traction vector for this stress tensor and surface normal is just

\begin{aligned}\mathbf{P} &=\frac{1}{{2}}\begin{bmatrix}1 & -2 & 2 \\ -2 & 3 & 1 \\ 2 & 1 & -1\end{bmatrix}\begin{bmatrix}\sqrt{2} \\ -1 \\ 1\end{bmatrix} \\ &=\frac{1}{{2}}\begin{bmatrix}\sqrt{2} + 2 + 2 \\ -2\sqrt{2} - 3 + 1 \\ 2\sqrt{2} - 1 -1\end{bmatrix} \\ &=\begin{bmatrix}\sqrt{2}/2 + 2 \\ -\sqrt{2} -1 \\ \sqrt{2} - 1\end{bmatrix}\end{aligned}

We also want a vector in the plane, and can pick

\begin{aligned}\boldsymbol{\tau} = \frac{1}{{\sqrt{2}}}\begin{bmatrix}0 \\ 1 \\ 1\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(4.25)

or

\begin{aligned}\boldsymbol{\tau}' = \begin{bmatrix}\frac{1}{{\sqrt{2}}} \\ \frac{1}{{2}} \\ -\frac{1}{{2}}\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(4.26)

It’s clear that either of these is normal to \hat{\mathbf{n}} (the first can also be computed by normalizing \hat{\mathbf{n}} \times \mathbf{e}_1, and the second with one round of Gram-Schmidt). However, neither of these vectors in the plane are particularly interesting since they are completely arbitrary. Let’s instead compute the projection and rejection of the traction vector with respect to the normal. We find for the projection

\begin{aligned}(\mathbf{P} \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}}&=\frac{1}{{4}}\left(\begin{bmatrix}\sqrt{2}/2 + 2 \\ -\sqrt{2} -1 \\ \sqrt{2} - 1\end{bmatrix}\cdot \begin{bmatrix}\sqrt{2} \\ -1 \\ 1\end{bmatrix} \right)\begin{bmatrix}\sqrt{2} \\ -1 \\ 1\end{bmatrix} \\ &=\frac{1}{{4}}\left( 1 + 2\sqrt{2}+\sqrt{2} +1 +\sqrt{2} - 1\right)\begin{bmatrix}\sqrt{2} \\ -1 \\ 1\end{bmatrix} \\ &=\frac{1}{{2}}\left( 1 + 4\sqrt{2}\right)\hat{\mathbf{n}}\end{aligned}

Our rejection, the component of the traction vector in the plane, is

\begin{aligned}(\mathbf{P} \wedge \hat{\mathbf{n}}) \hat{\mathbf{n}} &=\mathbf{P} - (\mathbf{P} \cdot \hat{\mathbf{n}})\hat{\mathbf{n}} \\ &=\frac{1}{{2}}\begin{bmatrix}\sqrt{2}/2 + 2 \\ -\sqrt{2} -1 \\ \sqrt{2} - 1\end{bmatrix}-\frac{1}{{4}}(1 + r \sqrt{2})\begin{bmatrix}\sqrt{2} \\ -1 \\ 1\end{bmatrix} \\ &=\frac{1}{{4}}\begin{bmatrix}\sqrt{2} \\ -3 \\ -5\end{bmatrix}\end{aligned}

This gives us a another vector perpendicular to the normal \hat{\mathbf{n}}

\begin{aligned}\hat{\boldsymbol{\tau}} = \frac{1}{{6}}\begin{bmatrix}\sqrt{2} \\ -3 \\ -5\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(4.27)

Wrapping up, we find the decomposition of the traction vector in the direction of the normal and its projection onto the plane to be

\begin{aligned}\mathbf{P} = \frac{1}{{2}}(1 + 4\sqrt{2}) \hat{\mathbf{n}}+\frac{3}{2} \hat{\boldsymbol{\tau}}.\end{aligned} \hspace{\stretch{1}}(4.28)

The components we can read off by inspection.

Problem Q4.

Statement

The stress tensor of a body is given by

\begin{aligned}\sigma =\begin{bmatrix}A \cos x & y^2 & C x \\ y^2 & B \sin y & z \\ C x & z & z^3\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(5.29)

Determine the constant A, B, and C if the body is in equilibrium.

Solution

In the absence of external forces our equilibrium condition was

\begin{aligned}\partial_k \sigma_{ik} = 0.\end{aligned} \hspace{\stretch{1}}(5.30)

In matrix form we wish to operate (to the left) with the gradient coordinate vector

\begin{aligned}0 &= \sigma \stackrel{ \leftarrow }{\boldsymbol{\nabla}} \\ &=\begin{bmatrix}A \cos x & y^2 & C x \\ y^2 & B \sin y & z \\ C x & z & z^3\end{bmatrix}\begin{bmatrix}\stackrel{ \leftarrow }{\partial}_x \\ \stackrel{ \leftarrow }{\partial}_y \\ \stackrel{ \leftarrow }{\partial}_z \\ \end{bmatrix} \\ &=\begin{bmatrix}\partial_x (A \cos x) + \partial_y(y^2) + \not{{\partial_z(C x)}} \\ \not{{\partial_x (y^2)}} + \partial_y(B \sin y) + \partial_z(z) \\ \partial_x (C x) + \not{{\partial_y(z)}} + \partial_z(z^3)\end{bmatrix} \\ &=\begin{bmatrix}-A \sin x + 2 y \\ B \cos y + 1 \\ C + 3 z^2 \end{bmatrix} \\ \end{aligned}

So, our conditions for equilibrium will be satisfied when we have

\begin{aligned}A &= \frac{2 y }{\sin x} \\ B &= -\frac{1}{\cos y} \\ C &= -3 z^2,\end{aligned} \hspace{\stretch{1}}(5.31)

provided y \ne 0, and y \ne \pi/2 + n\pi for integer n. If equilibrium is to hold along the y = 0 plane, then we must either also have A = 0 or also impose the restriction x = m \pi (for integer m).

A couple other mathematica notebooks

Some of the hand calculations done in this problem set I’ve confirmed using Mathematica. Those notebooks are available here

\begin{itemize}
\item

https://raw.github.com/peeterjoot/physicsplay/master/notes/phy454/continuumProblemSet1Q1.cdf

\item

https://raw.github.com/peeterjoot/physicsplay/master/notes/phy454/continuumProblemSet1Q2.cdf

\item

https://raw.github.com/peeterjoot/physicsplay/master/notes/phy454/continuumProblemSet1Q2animated.cdf

\item

https://raw.github.com/peeterjoot/physicsplay/master/notes/phy454/continuumProblemSet1Q3.cdf

\end{itemize}

These all require either a Mathematica client or the free Wolfram CDF player. Note that I haven’t figured out a way to get a browser based CDF player to play these without explicit download.

References

[1] L.D. Landau, EM Lifshitz, JB Sykes, WH Reid, and E.H. Dill. Theory of elasticity: Vol. 7 of course of theoretical physics. 1960.

Posted in Math and Physics Learning. | Tagged: , , , , , , , , , , , , , | 10 Comments »

« Previous Entries
 
Follow

Get every new post delivered to your Inbox.