An enumeration of other things I have written, including archives of all the individual pdfs that I have posted over the years along with my blog entries. All these pdfs are now stored directly on the new site. I will no longer be posting new content to any of my (three) google sites pages, nor will I be updating anything previously posted there. If you are looking for any corrections I may or may not have made, please look at the docs anchored off of the new blog.

A chronological listing of all the Mathematica notebooks I have written. The newest versions of these notebooks can still be found in my Mathematica github repository. A snapshot of each of these is now also available on the new site, so if you have the CDF plugin installed, these can now be examined by clicking on the links directly. Ironically, with chrome and my CDF installation, I’m able to view the .nb suffixed notebooks directly in the browser, but a click on any CDF (.cdf) notebook triggers a download?

Why a new blog?

Why after 611 blog posts on this wordpress.com hosted blog, dating all the way back to 2009, would I decide to ante-up and pay for hosting?

My primary motivation for this was truly geeky. I wanted the flexibility to be able to manage wordpress plugins (i.e. mathjax-latex and wolframcdf), and to also be able to put plain old html and arbitrary file content into the apache2 directory structure. I’ve wanted plain html hosting for a while, but made do with google sites (i.e. crappy but free). I’d also wanted to be able to use the wolfram CDF plugin on my blog, but also not enough to pay for it. However, once I tried mathjax-latex, I was sold. At least as a writer, compared to wp-latex, this “new way” completely kicks ass. It now takes much less time to produce posts with mathematics content, and requires far less scripting to convert from standalone latex. Unfortunately, one of the costs of this is pushed onto the reader, since it takes more time for mathjax-latex content to be formatted than the images produced by the standard wordpress latex plugin.

I tried out an amazon EC2 bitnami image for a while (amazon offers a free trial year to evaluate their offerings). That’s a flexible setup and offers direct access to the Linux VM, which is very nice. However, with an amazon EC2 image, I’m not really sure what I would end up paying. The charts seem somewhat vague, depending on future usage of both machine and storage. I would also have pay separately for a domain name, and pay separately for amazon hosting of the DNS entry.

I ended up deciding to use a go-daddy hosted wordpress instance, which is a flat rate service. It is less flexible than a godaddy standalone web-hosting environment, but also cheaper ($12 for the first year, including the domain name, and ~$50/year after that). It also looks like I can upgrade this to a more generic web hosting environment later if the cost of that seems justified. I’ll see first if only having sftp access to htdocs is enough of a major inconvenience to pay that additional yearly fee. EDIT: go-daddy either lied about their non-introductory rate, or increased their rates after I signed up. The renewal rates are about ~$100/year currently.

Configuring a custom MathJax configuration was a bit of a pain with only sftp access, mostly because I had to copy the MathJax tree, which was very slow for so many small files. I did that directory tree transfer with FileZilla since sftp ‘put –r’ appears to be busted. This MathJax setup was way easier on the EC2 since the ssh shell allowed for wget and local unzip directly from the apache2 htdocs tree. It’s a shame that the mathjax-latex plugin doesn’t allow the MathJax tree to be served from the default server (what the plugin settings calls the ‘MathJax CDN Service’). Logically, I’d like to be able to use that CDN service, but have my configuration file hosted locally. That config file (config/default.js) is a single small file, and is likely all that I’ll ever have to alter in that whole directory tree.

Solving for the flow below the stirring point.

Previously, we found the functional form for an azimuthal flow that has -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 above the bottom of the cup, and set the radius of that inner cylinder to as before, with the mixing occurring at an angular velocity of . We now want apply these boundary value constrains to the velocity function we’ve found for the steady state problem. Provided , this will have the form

where are the zeros of the order one Bessel function . Observe that our boundary value conditions of and are automatically satisfied. Note that because of the equality here, this form of solution is no good if we are mixing at the edge of the cup, so we require .

Our remaining boundary value condition , means that we have to solve

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

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

For a cup size of , stir radius of , stir depth , and angular velocity , we find

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 . 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 http://youtu.be/BS8XQdXljSk, 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.

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

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

Motivation and statement.

Viscous fluid is at rest in a two-dimensional channel between stationary rigid walls with . For a constant pressure gradient is imposed. Show that satifies

and give suitable initial and boundary conditions. Find in form of a Fourier series, and show that the flow approximates to steady channel flow when .

Solution

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

\begin{subequations}

\end{subequations}

Substitution of 2.2b into 2.2a gives us

Our equation to solve is therefore

This equation, for , allows for solutions

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

The no-slip conditions introduce boundary value conditions .

For we have

If we attempt separation of variables with , our equation takes the form

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

which has solution

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

Adding this we have , and subtracting gives us , 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

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

Applying the Navier-Stokes equation to this gives us

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

where our boundary value conditions are now given by

and

or
\begin{subequations}

\end{subequations}

Now we can apply separation of variables with , yielding

or

Here a positive constant has been used assuming that we want a solution that is damped with time.

Our solutions are

or

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

Rather remarkably, this Fourier series is actually a very good fit even after only a single term. Using the viscosity and density of water, , and (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).

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 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 http://youtu.be/0vZuv9HBtmo.

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 http://youtu.be/dDkx8iLwOew) 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.

Problem Q3 (revisited).

I’d produced the following sketches. For a higher viscosity bottom layer , 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 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

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 and respectively, and let the heights from the interface be and respectively.

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

Now it’s super easy to match the boundary conditions at and (the lower and upper walls respectively). Clearly the integration constants are just the velocities. Matching the tangential component of the traction vectors at we have

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 , fluid depths of and respectively, and viscosities and . We’ll write . We have a pair of Navier-Stokes equations to solve

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

Let’s require no components of the flows in the , or 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 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 we can restate Navier-Stokes in coordinate form as

\begin{subequations}

\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 we want to compute

So the tangential component of the traction vector is

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}

\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

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

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 our pressure is

so is constant with value

which fully determines the density of the upper surface

Matching the pressure between the two layers of fluids we have

so that our undetermined function

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

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 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

Observe that we have the usual 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 and

While more general solutions can likely be found, we will limit ourselves to looking only for solutions that are functions of . From our solution to the pressure part of the problem , we also see that the pressure term of 2.5b is killed. We are left with just

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

\begin{subequations}

\end{subequations}

The boundary conditions left to exploit are

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 we see immediately that we have . From the traction vector equality in the atmosphere, we have

or

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

and our traction vectors are

Matching both at the interface () gives us

We find

Let’s substitute this back for our first fluid velocity

As a check we see this is consistent with the previous calculation when . For our final integration constant we now find

So, finally, we have for the velocities

\begin{subequations}

\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.

Motivation.

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

Attempting to evaluate the integrals.

With

This is found to be

It is clear that we can’t evaluate this limit directly for since that gives us 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

and

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 coordinate integral first

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

where

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.