[cig-commits] r12357 - mc/2D/ConMan/trunk/doc
wei at geodynamics.org
wei at geodynamics.org
Mon Jun 30 10:59:49 PDT 2008
Author: wei
Date: 2008-06-30 10:59:49 -0700 (Mon, 30 Jun 2008)
New Revision: 12357
Added:
mc/2D/ConMan/trunk/doc/viscous-flow.tex
mc/2D/ConMan/trunk/doc/viscous-flow2b.tex
Log:
Supplemental material on 'Mantle Dynamics and Thermodynamics', and S. King's notes on Finite Element Methods for Viscous Flow.
Added: mc/2D/ConMan/trunk/doc/viscous-flow.tex
===================================================================
--- mc/2D/ConMan/trunk/doc/viscous-flow.tex (rev 0)
+++ mc/2D/ConMan/trunk/doc/viscous-flow.tex 2008-06-30 17:59:49 UTC (rev 12357)
@@ -0,0 +1,1624 @@
+\documentclass[12pt,twoside]{article}
+\usepackage{times}
+\begin{document}
+
+
+\title{EAS 591D Mantle Dynamics and Thermodyanmics}
+\author{Scott D. King}
+\date{Spring 2006}
+\maketitle
+
+\newcommand{\dx}{\delta x}
+\newcommand{\dz}{\delta z}
+\newcommand{\uxz}{u \left ( x , z \right )}
+\newcommand{\Uxz}{\vec u \left ( x , z \right )}
+\newcommand{\vxz}{v \left ( x , z \right )}
+\newcommand{\uxpz}{u \left ( x + \dx , z \right )}
+\newcommand{\vxzp}{v \left ( x , z + \dz \right )}
+\newcommand{\rxz}{\rho \left ( x , z \right )}
+\newcommand{\rxpz}{\rho \left ( x + \dx , z \right )}
+\newcommand{\rxzp}{\rho \left ( x , z + \dz \right )}
+\newcommand{\drdt}{{{\partial \rxz}\over{\partial t}}}
+\newcommand{\dudx}{{{\partial \uxz}\over{\partial x}}}
+\newcommand{\dvdz}{{{\partial \vxz}\over{\partial z}}}
+\newcommand{\drdx}{{{\partial \rxz}\over{\partial x}}}
+\newcommand{\drdz}{{{\partial \rxz}\over{\partial z}}}
+
+\section{Viscous Flow}
+\pagestyle{myheadings}
+\markboth{EAS 591d}{Mantle Dynamics and Thermodynamics}
+
+This course is designed to give the student an understanding of the
+basic equations for viscous flow and their solution on a computer. I assume
+that the student is reasonably familiar with vector calculus and linear algebra
+and is at least not afraid of partial differential equations. In addition, it
+will be very helpful to have a general familiarity and experience with Fortran
+programming. Even though the goal of this course is to have the student solving
+viscous flow problems on the computer, we will work out a number of analytic
+solutions. It is critical to have a working knowledge of simple analytic
+solutions, both to test one's computer code against and to build up one's
+intuition of viscous flow problems. There will be no exams for this class and
+most all of the assignments will have a computer aspect. The class will be very
+informal. Once we cover the basics, there may be a lot of discussion of the
+results of assignments in class. In some cases the student should be forewarned
+that not all of the answers to the assignments will be known before hand. The
+course will concentrate on the finite element method, since that is what I
+typically use in my research. However, if time and interests permit, we will
+look at some of the recent advances in high order finite differences and
+spectral methods.
+
+\begin{quotation}
+\noindent
+First Note: I will adopt an {\em Eulerian} frame of reference
+in this course. That means that we will consider our points of observation
+fixed and watch the fluid as it passes by them. This is in contrast to a {\em
+Lagrangian} framework (formulation), where instead of a fixed set of observation
+points, we track the path of a fixed set of fluid parcels. Both viewpoints have
+advantages and drawbacks. For most convection problems, tracking the paths of
+a set of fluid parcels quickly becomes an intractable problem. In some cases,
+we truly want to know the origin of a parcel of fluid. Later in the course we
+will discuss some methods for tracking parcels of fluid in an Eulerian system.
+\end{quotation}
+
+\begin{quotation}
+\noindent
+Second Note: Unless otherwise stated, I will be using a 2D
+Cartesian geometry. First, because it will make working problems on the
+blackboard easier and illustrates all of the basic principles (where 3D becomes
+important I will discuss it). Second, the basic form of ConMan is 2D Cartesian.
+\end{quotation}
+
+\subsection{Derive Equations}
+
+We begin with two basic physical principals; conservation of mass and
+conservation of momentum (this is truly a poor choice of name for the second
+principal, and we will see why in the following sections). There are a number
+of ways to derive the basic equations and several forms that we can write these
+equations in (e.g., stress and strain-rate or velocities). Because I come from
+a fluid dynamics background and because of the formulation used in ConMan, I
+will begin with the velocity formulation; however, I will do many of the
+derivations several ways because often different problems can be solved much
+more easily by using one form or the other. A familiarity with both forms is
+useful.
+
+This material is covered in a number of textbooks on fluid dynamics. There are
+several which I find quite useful:
+
+\begin{quotation}
+\noindent
+{\em Physical Fluid Dynamics,} D. J. Tritton, Oxford University Press, 1988.
+{\em The equations of motion are discussed in Chapter 5. Chapter 8 will also be
+quite useful for us. This is pretty much at the level of my notes.}
+\medskip
+
+\noindent
+{\em Fluid Mechanics,} L. D. Landau and E. M. Lifshitz, Pergamon Press, 1959.
+{\em The equations of motion are derived on the first few pages, actually
+Chapters 1 and 2 are quite useful.}
+\medskip
+
+\noindent
+{\em An Introduction to Fluid Dynamics,} G. K. Batchelor, Cambridge, 1967.
+{\em Chapter 3 is somewhat useful, but quickly gets into rotation. Batchelor
+is a classic fluids text.}
+\medskip
+
+\noindent
+{\em Geodynamics: Applications of Continuum Physics to Geological Problems,} D.
+L. Turcotte and G. Schubert, Wiley, 1982. {\em Chapter 6 covers fluid mechanics
+at a similar level to this class - they only consider viscous flow.}
+\medskip
+
+\noindent
+{\em Hydrodynamic and Hydromagnetic Stability,} S. Chandrasekhar, Oxford,
+1961. {\em This book is THE classic text for convection. The level of
+mathematical expertise required is high, but you should at least know of its
+existence.}
+\end{quotation}
+
+\subsubsection{Conservation of Mass}
+
+Let us begin by considering an infinitesimal box whose sides are of length $\dx$
+by $\dz$. We consider the mass flowing into and out of the box in a time
+$\delta t$.
+
+\begin{picture}(300,250)(-40,0)
+\put(75,25){\framebox(150,150){}}
+\put(45,100){\line(1,0){25}}
+\put(230,100){\line(1,0){25}}
+\put(150,0){\line(0,1){25}}
+\put(150,180){\line(0,1){25}}
+\put(35,110){$\uxz$}
+\put(235,110){$\uxpz$}
+\put(105,10){$\vxz$}
+\put(85,190){$\vxzp$}
+\put(10,85){\line(0,-1){60}}
+\put(5,100){$\dz$}
+\put(10,115){\line(0,1){60}}
+\put(135,240){\line(-1,0){60}}
+\put(150,235){$\dx$}
+\put(165,240){\line(1,0){60}}
+\end{picture}
+
+\begin{eqnarray}
+\drdt \dx \dz & = & \left [ \uxz \rxz - \uxpz \rxpz \right ]\dz \nonumber \\
+& & + \left [ \vxz \rxz - \vxzp \rxzp \right ]\dx \label{eq:cont1}
+\end{eqnarray}
+now I use the approximations
+\begin{eqnarray}
+\uxpz & = & \uxz + (\dx)\dudx \label{eq:uxpz} \\
+\vxzp & = & \vxz + (\dz)\dvdz \label{eq:vxzp}\\
+\rxpz & = & \rxz + (\dx)\drdx \label{eq:rxpz}\\
+\rxzp & = & \rxz + (\dz)\drdz \label{eq:rxzp}.
+\end{eqnarray}
+Using \ref{eq:uxpz}, \ref{eq:vxzp}, \ref{eq:rxpz}, and \ref{eq:rxzp}
+for the second and fourth terms on the right hand side of \ref{eq:cont1}, we get
+\begin{eqnarray}
+& \uxpz \rxpz & = \nonumber \\
+&\left (\uxz+(\dx)\dudx\right )\left (\rxz+(\dx)\drdx\right ) & = \nonumber
+\\
+& \uxz \rxz + \uxz \drdx \dx + \rxz \dudx \dx & + \nonumber \\
+& \dudx \drdx (\dx)^2
+\end{eqnarray}
+and
+\begin{eqnarray}
+& \vxzp \rxpz & = \nonumber \\
+&\left (\vxz+(\dz)\dvdz\right )\left (\rxz+(\dz)\drdz\right ) & = \nonumber
+\\
+& \vxz \rxz + \vxz \drdz \dz + \rxz \dvdz \dz & + \nonumber \\ &
+\dvdz \drdz (\dz)^2
+\end{eqnarray}
+Dropping terms of order $(\dx)^2$ and $(\dz)^2 $ and simplifying,
+equation~\ref{eq:cont1} becomes
+\begin{eqnarray}
+\drdt & = & - \left [ \uxz \drdx + \rxz \dudx \right ] \nonumber \\
+& & - \left [ \vxz \drdz + \rxz \dvdz \right ]
+\label{eq:cont2}
+\end{eqnarray}
+or in vector notation
+\begin{eqnarray}
+\drdt & = & - \rxz \nabla \cdot \Uxz - \Uxz \cdot \nabla \rxz \nonumber \\
+& = & - \nabla \cdot \left ( \rxz\Uxz \right )
+\label{eq:veccont2}
+\end{eqnarray}
+where
+\begin{equation}
+\Uxz = \left ( \uxz , \vxz \right ).
+\end{equation}
+
+Note a special case of interest is when the density $\rho$ is constant. In this
+case, equation~\ref{eq:veccont2} reduces to
+\begin{equation}
+\nabla \cdot \left ( \Uxz \right ) = 0.
+\label{eq:incompress-cont}
+\end{equation}
+
+A slightly cleaner version of the conservation of mass can be derived if we
+stick to vector notation throughout the derivation. Because I will present the
+momentum equation fully in vector notation, it is worth going through the
+derivation again:
+
+\begin{picture}(300,100)(-40,-20)
+\put(150,25){\circle{40}}
+\put(155,30){\line(1,1){20}}
+\put(155,35){$\vec u$}
+\put(150,40){\line(0,1){20}}
+\put(130,45){$\vec dS$}
+\end{picture}
+
+We will consider some arbitrary volume of fluid, $V$, with a surface, $S$.
+We will define $d\vec S$, a vector with an outward pointing normal to the surface
+$S$. Then, the rate of loss of mass from $V$ will equal
+\begin{equation}
+\oint \rho \vec u \cdot d\vec S.
+\label{eq:rml}
+\end{equation}
+Transforming equation~\ref{eq:rml} into an integral over volume,
+\begin{equation}
+\int \nabla \cdot (\rho \vec u ) dV.
+\label{eq:rml2}
+\end{equation}
+The total mass in volume $V$ is given by
+\begin{equation}
+\int \rho dV,
+\end{equation}
+hence,
+\begin{equation}
+{d\over dt}\int \rho dV = -\int \nabla \cdot (\rho \vec u ) dV.
+\end{equation}
+Because we consider a constant volume
+\begin{equation}
+\int {{\partial \rho}\over{\partial t}} dV = -\int \nabla \cdot (\rho \vec
+u ) dV.
+\label{eq:cont3}
+\end{equation}
+The only way equation~\ref{eq:cont3} can be valid over the volume $V$ is if at
+every point
+\begin{equation}
+{{\partial \rho}\over{\partial t}} + \nabla \cdot (\rho \vec u ) = 0.
+\label{eq:cont4}
+\end{equation}
+Which should look familiar.
+
+\subsubsection{Conservation of Momentum}
+
+The next conservation equation, often referred to as the conservation of
+momentum, is not so easy to derive. I am also going to drop the $(x,z)$
+notation from my variables, but remember that in general velocities and
+densities can be a function of their location in space and time.
+
+We will consider some arbitrary volume of fluid, $V$, with a surface, $S$,
+as in the picture above. The total force, $T_f$ acting on this volume of fluid
+is equal to the integral of the pressure, $p$, taken over the surface of the
+fluid.
+\begin{equation}
+- \oint p dS = T_f
+\label{eq:pres1}
+\end{equation}
+Transforming \ref{eq:pres1} to a volume integral
+\begin{equation}
+- \oint p dS = - \int \nabla p dV.
+\label{eq:pres2}
+\end{equation}
+
+At this point, I am going to neglect possible effects due to dissipation of
+energy of the fluid. We will build that in at the next step. From
+equation~\ref{eq:pres2}, we see that the fluid surrounding any volume element
+$dV$ exerts a force on that element of $- \nabla p$. That force must be balanced by a change in acceleration
+of the fluid according to Newton's second law.
+\begin{equation}
+{{D (\rho\vec u)} \over {D t}} = - \nabla p
+\label{eq:euler1}
+\end{equation}
+Notice that
+\begin{equation}
+{{D (\rho\vec u)}\over{D t}}=\rho{{D\vec u}\over{Dt}} + \vec u {{D \rho}
+\over {D t}}
+\label{eq:totalmomentum}
+\end{equation}
+however the second term in Equation~\ref{eq:totalmomentum} is zero from
+Equation~\ref{eq:cont4} The only physical mechanisms for this would be changes
+in phase or if velocities approached the speed of light and both of these cases
+would (and can) be handled as special cases. However, in most problems of
+geophysical interest, the fluids are not moving at anything approaching the
+speed of light. A fluid with a change in phase could possibly be a special
+topic to be covered later.
+
+Note that ${{D \vec u} \over {D t}}$ is not the rate of change of the fluid
+velocity at a fixed point in space, but the rate of change of a given fluid
+particle as it moves about in space (that is why I did not write it in the
+usual way). However, because we are in an Eulerian framework, it has to be
+expressed in terms of quantities referring to a fixed point in space. To
+accomplish this, consider an arbitrarily small parcel of fluid moving between two
+points (arbitrarily close), denoted by the vector $d \vec r$. The parcel has two
+components to its change in velocity $d\vec u$, namely the change due to the
+change during the time $\delta t$ at a point fixed in space, and the
+difference between the velocities, at the same instant in time, at the two
+points $d \vec r$ apart. We can write this as
+\begin{equation}
+d \vec u = \left ({{\partial \vec u} \over {\partial t}} \right ) \delta t +
+(d\vec r \cdot \nabla) \vec u,
+\end{equation}
+or dividing both sides by $\delta t$,
+\begin{equation}
+{{D \vec u} \over {D t}} = {{\partial \vec u} \over {\partial t}} +
+(\vec u \cdot \nabla) \vec u.
+\label{eq:total-deriv}
+\end{equation}
+Substituting \ref{eq:total-deriv} into \ref{eq:euler1}, we get
+\begin{equation}
+{{\partial \vec u} \over {\partial t}} + (\vec u \cdot \nabla) \vec u =
+- {1\over\rho}\nabla p.
+\label{eq:euler-final}
+\end{equation}
+Equation~\ref{eq:euler-final} is the equation of motion of a fluid, with no
+dissipation. It was first derived by L. Euler in 1755 and is known as {\em
+Euler's equation}. If the fluid is in a gravitational field we need to add an
+additional body force to the right hand side $\rho \vec g$ where $\vec g$ is
+the acceleration due to gravity. This modifies \ref{eq:euler-final},
+\begin{equation}
+{{\partial \vec u} \over {\partial t}} + ( \vec u \cdot \nabla ) \vec u =
+- {1\over\rho}\nabla p + \vec g.
+\label{eq:euler-final2}
+\end{equation}
+This is the basic equation governing fluid motion. It is the fundamental
+equation for atmospheric science, hydrology, and most other fluid sciences,
+even though the way some of these fields write their basic equations many
+appear in a different form than \ref{eq:euler-final2}.
+
+This is a good place to mention the {\em Boussinesq Approximation,} or what is
+sometimes called {\em incompressible flow.} In the Boussinesq approximation, we
+assume that the density is constant everywhere, except allowing the density to
+vary in the buoyancy term of equation~\ref{eq:euler-final2}. Thus, the
+continuity equation reduces to
+\ref{eq:incompress-cont} and Euler's equation becomes
+
+\begin{equation}
+{{\partial \vec u} \over {\partial t}} + ( \vec u \cdot
+\nabla ) \vec u = - {1\over\rho_o}\nabla p + {\rho\over\rho_o}\vec g.
+\label{eq:euler-final3}
+\end{equation}
+You might wonder about the validity of this approximation. However, you can
+convince yourself that if $||{{D \vec u} \over {D t}}||$ is very small compared
+to $||\vec g||$, as it most certainly is in most geologically relevant cases,
+then the error made by neglecting changes in density in this term will be
+small. The term ${\rho\over\rho_o}\vec g$ in equation~\ref{eq:euler-final3} is
+the connection between density perturbations from temperature, which drive flow
+in the thermal convection problem. Without variations in density, the flow is
+only driven by pressure variations and/or velocity boundary conditions.
+
+ConMan, and many other fluid flow codes use the Boussinesq approximation;
+however, a number of researchers have shown that in mantle convection
+compressibility effects (changes in density) can have important effects on the
+flow, so compressible flow is becoming more and more common. If time permits,
+we will examine compressibility in greater detail later in the course.
+
+\subsubsection{Conservation of Momentum: Viscous Dissipation}
+
+We now need to add the viscous forces acting on the element. For this we need
+to recall some basic principals of stress and strain-rate. If you need to
+refresh your memory Turcotte and Schubert Chapter 2 is a pretty good place to
+look. If you look in Tritton or Landau and Lifshitz, you will be pretty
+disappointed, probably because most fluid dynamics problems are not dominated
+by stress and strain-rate.
+
+It is probably worth a brief digression to remind everyone of the fundamentals
+of stress and strain (strain-rate). If a solid is {\em elastic,} then
+stress is linearly related to strain. If a material is a {\em Newtonian
+fluid,} then stress is linearly related to strain-rate. I find it easiest to
+think in terms of elastic solids and then take the time derivative to get
+strain-rates. I can't think of any particular reason why this is not a valid
+thing to do.
+
+Lets consider a 2D box, whose dimensions are $\dx$ and $\dz$. Lets squeeze the
+box so that it's new dimensions are $\dx^\prime$ $\dz^\prime$, we could imagine
+doing this by increasing hydrostatic pressure. Notice that in general the box
+may also translate, i.e., the center of mass moves. This translation, should
+not effect the strain in any way. In addition, there could be rigid body
+rotation, this also should not effect the strain.
+
+\begin{picture}(300,250)(-40,0)
+\put(75,25){\framebox(150,150){}}
+\put(115,55){\framebox(130,130){}}
+
+\put(10,85){\line(0,-1){60}}
+\put(5,100){$\dz$}
+\put(10,115){\line(0,1){60}}
+\put(135,240){\line(-1,0){60}}
+\put(150,235){$\dx$}
+\put(165,240){\line(1,0){60}}
+\put(30,105){\line(0,-1){50}}
+\put(25,120){$\dz^\prime$}
+\put(30,135){\line(0,1){50}}
+\put(165,220){\line(-1,0){50}}
+\put(180,215){$\dx^\prime$}
+\put(195,220){\line(1,0){50}}
+\end{picture}
+
+Now I should be careful in my use of notation, because we are doing an elastic
+problem, so velocities become displacements, to remind you of that, I'm going to
+use $w$ for displacement; however, lets charge merrily along... Lets assume
+the lower left corner of the box is at $(x,z)$ and after the ``strain event''
+is at
+$(x^\prime,z^\prime)$. Then the {\em displacement} of the lower left corner
+will be
+\begin{eqnarray}
+w_x (x,z) = x - x^\prime \\
+w_z (x,z) = z - z^\prime
+\end{eqnarray}
+and for the lower right hand corner will be
+\begin{equation}
+w_x (x + \dx,z) = (x + \dx) - (x^\prime + \dx^\prime )
+\end{equation}
+and for the upper left hand corner
+\begin{equation}
+w_z (x ,z + \dz) = (z + \dz) - (z^\prime + \dz^\prime ).
+\end{equation}
+Now, we once again do our usual trick, using a first-order Taylor expansion for
+the terms $w_x (x + \dx,z)$ and $w_z (x ,z + \dz)$, since $\dx$ and $\dz$ are
+infinitesimal.
+\begin{eqnarray}
+w_x (x + \dx,z) & = & w_x (x,z) + {{\partial w_x }\over{\partial x}} \dx \\
+w_z (x ,z + \dz) & = & w_x (x,z) + {{\partial w_z }\over{\partial z}} \dz
+\end{eqnarray}
+by substituting and rearranging, we can get
+\begin{eqnarray}
+\dx & = & \dx^\prime + {{\partial w_x }\over{\partial x}} \dx \\
+\dz & = & \dz^\prime + {{\partial w_z }\over{\partial z}} \dz.
+\end{eqnarray}
+Thus, if we define strain as the change in length divided by the original
+length, we get
+\begin{eqnarray}
+\epsilon_{xx} & = & {{\dx -\dx^\prime}\over \dx} = {{\partial w_x
+}\over{\partial x}}
+\\
+\epsilon_{zz} & = & {{\dz -\dz^\prime}\over \dz} = {{\partial w_z }\over{\partial
+z}}.
+\end{eqnarray}
+
+This simple little analysis only dealt with normal strains, that is shortening
+(or lengthening) in the direction parallel to the strain. We now need to look
+at shear strains. For that we need the figure below.
+
+\begin{picture}(300,250)(-40,0)
+\put(25,25){\framebox(150,150){}}
+\put(25,25){\line(4,1){150}}
+\put(25,25){\line(1,4){37.5}}
+\put(175,62.5){\line(1,4){37.5}}
+\put(62.5,175){\line(4,1){150}}
+\put(10,85){\line(0,-1){60}}
+\put(5,100){$\dz$}
+\put(10,115){\line(0,1){60}}
+\put(85,240){\line(-1,0){60}}
+\put(100,235){$\dx$}
+\put(115,240){\line(1,0){60}}
+\end{picture}
+
+I will define the shear strain as one-half the decrease in the right angle in
+the lower left corner. If I define $\phi_1$ and $\phi_2$ as the angles
+through which the sides of the original quadrilateral are rotated,
+then
+\begin{equation}
+\epsilon_{xz} = -{1\over2} ( \phi_1 + \phi_2 )
+\end{equation}
+Notice that by this definition, the shear strain is symmetric, i.e.,
+\begin{equation}
+\epsilon_{xz} = \epsilon_{zx}.
+\end{equation}
+Further, we can relate $\phi_1$ and $\phi_2$ to the displacements and use the
+fact that the displacements are infinitesimal, so that we can use $\tan \alpha
+= \alpha $
+\begin{eqnarray}
+\tan ~ \phi_1 & = & {{-w_z (x + \dx, z)}\over \dx} = \phi_1 \\
+\tan ~ \phi_2 & = & {{-w_x (x , z + \dz )}\over \dz} = \phi_2
+\end{eqnarray}
+Once again, we will use our first-order Taylor expansion trick again.
+Jumping straight to the result, the relation between shear strain and the
+derivatives of displacement
+\begin{equation}
+\epsilon_{xz} = {1\over 2} \left ( {{\partial w_z}\over {\partial x}} +
+{{\partial w_x}\over {\partial z}} \right ).
+\end{equation}
+
+For completeness, shear strain can lead to a solid body rotation. To see this,
+we need the definition of the solid body rotation about the $y$ axis
+\begin{equation}
+\omega_{y} = -{1\over2} ( \phi_1 - \phi_2 ) = {1\over 2} \left ( {{w_z}\over
+\dx} - {{w_x}\over \dz} \right )
+\end{equation}
+Using the relations above, you can write
+\begin{eqnarray}
+\phi_1 = - (\epsilon_{xz} + \omega_y ) \\
+\phi_2 = - (\epsilon_{xz} - \omega_y )
+\end{eqnarray}
+For the case $\phi_1 = \phi_2$, we have pure shear, and $\omega_y = 0$. For the
+case where one of the angles is zero, we have both a shear strain and a
+rotation.
+
+To convince yourself that you understand what is going on, it would be
+wise to look at a few simple problems.
+\begin{quote}
+{\bf Problem 1} Add two simple shear problems together, one in the $x$
+direction and one in the $z$ direction, both with the same angle $\alpha$.
+Show the resulting shear stress and rotation.
+\end{quote}
+
+\begin{quote}
+{\bf Problem 2} Show for the triangle below that
+\begin{equation}
+\tau_{x^\prime x^\prime} = \tau_{xx} \cos^2 \theta + \tau_{zz} \sin^2 \theta +
+\tau_{xz} \sin 2\theta
+\end{equation}
+(Note, I drew it like $\theta = \pi/4$ but you can not assume this!)
+\end{quote}
+
+\begin{picture}(300,180)(-40,0)
+\put(125,25){\line(0,1){150}}
+\put(125,25){\line(1,0){150}}
+\put(125,175){\line(1,-1){150}}
+\put(75,100){\line(1,0){20}}
+\put(80,105){$\tau_{xx}$}
+\put(200,0){\line(0,1){20}}
+\put(190,10){$\tau_{zz}$}
+\put(200,100){\line(1,1){20}}
+\put(210,110){$\tau_{x^\prime x^\prime}$}
+\put(105,80){\line(0,1){40}}
+\put(110,85){$\tau_{xz}$}
+\end{picture}
+
+\vfill\eject
+
+And now that our brief tour of strain is over, we will go back to fluids. In
+everything that follows, I will try and stick to the following convention,
+deviatoric stress will be denoted by $\tau$, deviatoric strain-rate by $\dot
+\epsilon$ and total stress by $\sigma$.
+
+\newcommand{\txxxz}{\tau_{xx} \left ( x , z \right )}
+\newcommand{\txxxpz}{\tau_{xx} \left ( x +\delta x , z \right )}
+\newcommand{\tzzxz}{\tau_{zz} \left ( x , z \right )}
+\newcommand{\tzzxpz}{\tau_{zz} \left ( x + \dx , z \right )}
+\newcommand{\tzzxzp}{\tau_{zz} \left ( x , z + \dz \right )}
+\newcommand{\txzxz}{\tau_{xz} \left ( x , z \right )}
+\newcommand{\txzxpz}{\tau_{xz} \left ( x + dx , z \right )}
+\newcommand{\tzxxz}{\tau_{zx} \left ( x , z \right )}
+\newcommand{\tzxxzp}{\tau_{zx} \left ( x , z + \dz \right )}
+
+\begin{picture}(300,250)(-40,-20)
+\put(75,25){\framebox(150,150){}}
+\put(40,100){\line(1,0){25}}
+\put(235,100){\line(1,0){25}}
+\put(150,-10){\line(0,1){25}}
+\put(150,190){\line(0,1){25}}
+\put(30,110){$\txxxz$}
+\put(240,110){$\txxxpz$}
+\put(95,10){$\tzzxz$}
+\put(75,190){$\tzzxzp$}
+
+\put(10,85){\line(0,-1){60}}
+\put(5,100){$\dz$}
+\put(10,115){\line(0,1){60}}
+\put(135,240){\line(-1,0){60}}
+\put(150,235){$\dx$}
+\put(165,240){\line(1,0){60}}
+\end{picture}
+
+The first thing to notice is that if there is going to be no net torque,
+then
+\begin{equation}
+\txzxz = \tzxxz.
+\end{equation}
+Next, following what we did in the case of conservation of mass,
+the net viscous force in the x direction per unit area is
+\begin{equation}
+{{\left [ \txxxpz - \txxxz \right ] \dz}\over \dx\dz} +
+{{\left [ \tzxxzp - \tzxxz \right ] \dx}\over \dx\dz} = 0.
+\label{eq:xxstress}
+\end{equation}
+Similarly for the $z$ direction, we have
+\begin{equation}
+{{\left [ \tzzxpz - \tzzxz \right ] \dx}\over \dx\dz} +
+{{\left [ \txzxpz - \txzxz \right ] \dz}\over \dx\dz} = 0.
+\label{eq:zzstress}
+\end{equation}
+
+Once again we do our Taylor expansion, for the $\txxxpz$ type terms and
+equation~\ref{eq:xxstress} reduces to
+\begin{equation}
+{{\partial\txxxz}\over{\partial x}} + {{\partial\tzxxz}\over{\partial z}} =0
+\label{eq:xxstress2}
+\end{equation}
+and equation~\ref{eq:zzstress} reduces to
+\begin{equation}
+{{\partial\tzzxz}\over{\partial z}} + {{\partial\txzxz}\over{\partial x}} =0.
+\label{eq:zzstress2}
+\end{equation}
+
+We can write equations~\ref{eq:xxstress2} and~\ref{eq:zzstress2} in a general
+tensor notation as
+\begin{equation}
+\tau_{ij,i} = 0
+\end{equation}
+where repeated subscripts mean to sum the terms over, in this case $i=1,2$ and
+$j=1,2$.
+
+For an ideal {\em Newtonian} viscous fluid, viscous stresses are linearly
+proportional to velocity gradients (by definition). This leads us to
+\begin{eqnarray}
+\tau_{xx} & = & 2 \mu {{\partial u}\over{\partial x}} \\
+\tau_{zz} & = & 2 \mu {{\partial v}\over{\partial z}} \\
+\tau_{xz} & = & \tau_{zx} = \mu \left ( {{\partial u}\over{\partial z}} +
+{{\partial v}\over{\partial x}} \right )
+\end{eqnarray}
+where $\mu$ is the dynamic viscosity (Pa s). Recall that a Pascal is ${\rm
+kg} / {\rm m s^2}$ is more fundamental units. Notice that in
+$ij$ notation we could write
+\begin{equation}
+\tau_{ij} = \mu \dot\epsilon_{ij}
+\end{equation}
+Substituting for velocity derivatives in equation~\ref{eq:xxstress2}, we get
+\begin{equation}
+{ {\partial \left (2 \mu {{\partial u}\over{\partial x}} \right ) } \over
+{\partial x} } + { {\partial \left ( \mu \left ( {{\partial u}\over{\partial z}}
++ {{\partial v}\over{\partial x}} \right ) \right )} \over {\partial z} } = 0
+\end{equation}
+and a similar messy equation for~\ref{eq:zzstress2}. Notice that {\em if and
+only if} the viscosity is a constant, we can simplify the above to
+\begin{equation}
+2 \mu {{\partial^2 u}\over{\partial x^2}} +
+\mu \left ( {{\partial^2 u}\over{\partial z^2 }} +
+{{\partial^2 v}\over{\partial x \partial z}} \right ) = 0.
+\end{equation}
+Now we return to the conservation of mass equation~\ref{eq:incompress-cont} for
+the incompressible case, and differentiate each side by $\partial x$
+\begin{equation}
+{{\partial^2 v}\over{\partial x \partial z}} = - {{\partial^2 u}\over{\partial
+x^2 }}
+\end{equation}
+so in this case, we arrive at
+\begin{equation}
+\mu \left ( {{\partial^2 u}\over{\partial x^2 }} + {{\partial^2
+u}\over{\partial z^2}} \right ) = 0.
+\end{equation}
+In a exactly the same way, we can get the term for the viscous stress in the
+$z$ direction
+\begin{equation}
+\mu \left ( {{\partial^2 v}\over{\partial x^2 }} + {{\partial^2
+v}\over{\partial z^2}} \right ) = 0.
+\end{equation}
+or in vector notation
+\begin{equation}
+\mu \nabla^2 \vec u = 0
+\end{equation}
+Notice in vector notation, we can write the case where $\mu$ is not constant
+more cleanly,
+\begin{equation}
+\nabla ( \mu \nabla \vec u) = 0.
+\end{equation}
+
+If we add this term to our Euler equation, from back on
+equation~\ref{eq:euler-final2}, we have
+\begin{equation}
+{{\partial \vec u} \over {\partial t}} + ( \vec u \cdot
+\nabla ) \vec u = - {1\over\rho}\nabla p + \vec g + {1\over\rho}
+\nabla ( \mu \nabla \vec u).
+\label{eq:navier-stokes1}
+\end{equation}
+
+Now we have the most general equation for describing a the motion of an
+irrotational, non-magnetic fluid. If you want to see the effects of rotation
+and magnetic fields, go look at Chandrasekhar. Sometimes the term $\mu/\rho$ is
+called the dynamic viscosity, with units of ${\rm m^2}/ s$. Thus, special care
+needs to be taken when discussing viscosity.
+
+\vfill\eject
+
+\subsection{Non-Dimensionalization}
+
+What is non-dimensionalization? We rewrite the equations in terms of
+dimensionless variables by dividing out a constant characteristic size for the
+length, time, velocity, etc in the equations. Because these are constant, we
+can take them outside the derivatives and group them into parameters. These
+parameters allow us to see the relationship between the variables and to scale
+problems from the laboratory or computer, to Earth scales.
+
+Why do we want to non-dimensionalize equations? Often,
+the equations we want to solve in fluid problems can not be solved with all the
+terms present, sometimes terms that make an equation intractable do not
+contribute to the physics of the problem we want to solve. Also, by using
+dimensional analysis, we can often gain insight into what physical processes
+control (dominate) the problems we are interested in. As we will see later,
+this is also something important to do on a computer where we have a finite
+precision.
+
+It is helpful to think of our variables in terms of the fundamental units of
+mass, $M$, length, $L$, and time, $T$
+\begin{eqnarray*}
+x & \rightarrow & L \\
+t & \rightarrow & T \\
+\rho & \rightarrow & {M\over L^3} \\
+\rho \vec g & \rightarrow & {{M}\over{L^2 T^2}} \\
+\tau & \rightarrow & {M\over {LT^2}} \\
+p & \rightarrow & {M\over {LT^2}} \\
+\dot \epsilon & \rightarrow & {1\over T} \\
+\mu & \rightarrow & {{M}\over{LT}}
+\end{eqnarray*}
+
+Lets begin in the following way, I will use primed variables to
+represent the non-dimensional form of the variable. Lets think of solving for
+the motion of the fluid in a box whose depth is $D$. We need a characteristic
+time scale for time. To do this, our first assumption will be that we know
+something about the characteristic velocity of the flow.
+
+This leads to the following relations between the dimensional and
+non-dimensional variables
+\begin{eqnarray}
+x & = & x^\prime D \\
+u & = & u^\prime {V} \\
+t & = & t^\prime {D\over V} \\
+p & = & P_o + p^\prime {\mu V\over D}
+\end{eqnarray}
+Here I have pulled out a background pressure, $P_o$, because it is only
+gradients in pressure that drive flow (see~\ref{eq:navier-stokes1}). Keep in
+mind that derivatives must be dimensionless too, so
+\begin{equation}
+\nabla = \nabla^\prime {1\over D}
+\end{equation}
+Plugging these into equation~\ref{eq:navier-stokes1}, we get
+\begin{equation}
+{{\partial \vec u^\prime} \over {\partial t^\prime}}{{V}\over{D\over V}}
+ + ( \vec u^\prime \cdot \nabla^\prime ) \vec u^\prime {V^2\over D}
+= -{1\over\rho}\nabla^\prime p^\prime {1\over D} {\mu V\over D}
++ {{\delta\rho}\over\rho}\vec g^\prime {{\mu V}\over D^2} +
+{\mu\over\rho} \left (\nabla^\prime \right )^2 \vec u^\prime {V\over D^2}.
+\label{eq:navier-stokes2}
+\end{equation}
+Dividing through by ${{V^2}\over{D}}$
+\begin{equation}
+{{\partial \vec u^\prime} \over {\partial t^\prime}}
+ + ( \vec u^\prime \cdot \nabla^\prime ) \vec u^\prime =
+-\nabla^\prime p^\prime {\mu \over {\rho V D}}
++ \delta\rho \vec g^\prime {\mu \over {\rho V D}}
++ {\mu\over\rho} \left (\nabla^\prime \right )^2 \vec u^\prime {1\over {V D}}.
+\label{eq:navier-stokes-nd}
+\end{equation}
+which we can rewrite as
+\begin{equation}
+{{\partial \vec u^\prime} \over {\partial t^\prime}}
+ + ( \vec u^\prime \cdot \nabla^\prime ) \vec u^\prime
+= -\nabla^\prime p^\prime {1\over {\rm Re}}
++ \delta\rho \vec g^\prime {1\over {\rm Re}}
++ {1\over {\rm Re}} \left (\nabla^\prime \right )^2 \vec u^\prime.
+\label{eq:navier-stokes-nd2}
+\end{equation}
+where
+\begin{equation}
+Re = {{\rho V D}\over\mu}
+\end{equation}
+is the Reynolds number. There are several interesting things to think about
+with the Reynolds number. First examining the terms, you see that the
+Reynolds number is dimensionless. In fluid dynamics, there are a set,
+actually a whole bunch, of dimensionless numbers that govern the dynamics of
+fluid motions. Next, notice that any combination of $\rho, V, D, {\rm and}
+\mu$ which give the same Reynolds number will give the same answer to the
+equations. Hence, if we solve the dimensionless problem for a given Reynolds
+number, we have a solution for a whole class of problems where $\rho, V, D, {\rm
+and} \mu$ give the same Reynolds number.
+
+Lets look at some specific examples of flows:
+
+\begin{itemize}
+\item Motion of the atmosphere:
+\begin{eqnarray*}
+\rho & = & 1.3~{{\rm kg}\over{\rm m^3}} \\
+V & = & 20.0~{{\rm km}\over{\rm hour}} \\
+D & = & 10~{\rm km} \\
+\mu & = & 1.78 \time 10^{-5}~{\rm Pa~s} \\
+Re & = & {{1.3 \times 2.0 \times 10^4 \times 2.7 \times 10^{-4} \times 10^4 }
+\over {1.78 \times 10^{-5}} } \\
+ & = & 6 \times 10^{9}
+\end{eqnarray*}
+
+\item Water in a stream:
+\begin{eqnarray*}
+\rho & = & 10^3~{{\rm kg}\over{\rm m^3}} \\
+V & = & 1.0~{{\rm m}\over{\rm s}} \\
+D & = & 10~{\rm m} \\
+\mu & = & 1.14\times 10^{-3}~{\rm Pa~s} \\
+Re & = & {{10^3 \times 1.0 \times 10.} \over {1.14 \times 10^{-3}} } \\
+ & = & 10^7
+\end{eqnarray*}
+
+\item A surging glacier:
+\begin{eqnarray*}
+\rho & = & 900~{{\rm kg}\over{\rm m^3}} \\
+V & = & 10.0~{{\rm m}\over{\rm day}} \\
+D & = & 100.0~{\rm m} \\
+\mu & = & 1.0\times 10^{6}~{\rm Pa~s} \\
+Re & = & {{900 \times 10.0 \times 1 \times 10^{-5} \times 100. } \over {1.0
+\times 10^{10} }} \\
+ & = & 9 \times 10^{-10}
+\end{eqnarray*}
+
+\item Motion of the Mantle Resulting from Plate Motions:
+\begin{eqnarray*}
+\rho & = & 3300~{{\rm kg}\over{\rm m^3}} \\
+V & = & 10.0~{{\rm cm}\over{\rm yr}} \\
+D & = & 3000.0~{\rm km} \\
+\mu & = & 1.0\times 10^{21}~{\rm Pa~s} \\
+Re & = & {{3300 \times 10.0 \times 3 \times 10^{-10} \times 3\times 10^6 }
+\over {1.0 \times 10^{21} }} \\
+ & = & 3 \times 10^{-20}
+\end{eqnarray*}
+
+\end{itemize}
+
+Notice that when $Re \rightarrow 0$ equation~\ref{eq:navier-stokes-nd2} reduces
+to {\em Stokes equation}
+\begin{equation}
+0 = -\nabla^\prime p^\prime + \delta\rho \vec g^\prime + \left
+(\nabla^\prime \right )^2 \vec u^\prime.
+\label{eq:stokes-nd}
+\end{equation}
+
+Considering the other end-member, when $Re \rightarrow \infty ,$
+equation~\ref{eq:navier-stokes-nd2} almost reduces to Euler's equation. If we
+have chosen to scale the pressure by $\rho V^2$ the Reynolds number did not show
+up in front of the pressure term, then equation~\ref{eq:navier-stokes-nd2} would
+have reduced to Euler's equation. This scaling ambiguity, suggests that the
+pressure gradient term should always be left in the equations unless one has a
+good physical mechanism to get rid of it.
+
+Luckily for Earth Scientists, many of the problems we want to solve are like
+the problems I illustrated above--the Reynolds number is either very larger
+or very small. In either case, one group of terms can be
+neglected, greatly simplifying the problem. Note that we did not consider the
+effect of rotation in our derivation of the equations. The Earth is a
+rotating body and all of the fluids moving around on and in Earth are subject
+to the Coriolis force. I don't want to go back and add that contribution, but
+I will state that when you go through the non-dimensionalization, you get a
+Reynolds number in front of that term as well. Because the focus of this
+course is on Stokes equation, or viscous flow, we will not mention the
+effect of rotation any further.
+
+There are several things to think about. Stokes equation assumes
+viscous stresses are balanced by buoyancy and pressure gradients. There is no
+inertia in the problem. If you change the force, the fluid changes
+instantaneously in response to the change in force. There is also no time
+dependence in Stokes equation. Given a sent of boundary conditions for
+velocity, Stokes equation gives a steady solution (time enters no where in the
+equation). Time can only enter in through the changing of the buoyancy terms
+$\delta \rho$ or time-dependent boundary conditions.
+
+At this point, you should be wondering what happens when we consider a case
+where we don't have any idea what characteristic velocity we should use in the
+scaling. That is the case, for example, in convection problems. We will
+consider that case next class.
+
+Now lets look at another way to scale time. For this problem, we again
+consider the motion of the fluid in a box whose depth is $D$. Now, we
+will assume that the fluid is convecting. That means, the $\delta \rho$ term
+in the buoyancy term is the result of temperature gradients in the fluid. The
+temperatures will be advected by the flow, making the problem non-linear.
+This adds a number of complications, which for now we will not worry about.
+However, we don't have any idea what the velocity should be, so we can't use a
+characteristic velocity as in the Reynolds number case above. Another time
+scale of interest will be related to the diffusion of heat. The parameter of
+interest is the thermal diffusivity, $\kappa$. By inspecting the units of
+thermal diffusivity (${\rm m^2 / sec}$, we can get a time scale of $D^2/\kappa$.
+
+This leads to the following relations between the dimensional and
+non-dimensional variables
+\begin{eqnarray}
+x & = & x^\prime D \\
+t & = & t^\prime {D^2 \over \kappa} \\
+u & = & u^\prime {\kappa \over D} \\
+p & = & P_o + p^\prime {{\mu \kappa}\over D^2}
+\end{eqnarray}
+Once again, I have pulled out a background pressure, $P_o$, because it
+is only gradients in pressure that drive flow (see~\ref{eq:navier-stokes1}).
+
+Plugging these into equation~\ref{eq:navier-stokes1}, we get
+\begin{equation}
+{{\partial \vec u^\prime} \over {\partial t^\prime}}
+{{\kappa\over D}\over{D^2\over \kappa}}
++ (\vec u^\prime \cdot \nabla^\prime ) \vec u^\prime {{\kappa^2}\over D^3} =
+-{1\over\rho}\nabla^\prime p^\prime {1\over D} {{\mu\kappa}\over D^2} +
+{{\delta\rho}\over\rho}\vec g^\prime {{\mu \kappa}\over D^3} + {\mu\over\rho}
+\left (\nabla^\prime \right )^2 \vec u^\prime {\kappa\over D^3}.
+\label{eq:navier-stokes4}
+\end{equation}
+Dividing Equation~\ref{eq:navier-stokes4} by $\kappa^2\over D^3$ yields,
+\begin{equation}
+{{\partial \vec u^\prime} \over {\partial t^\prime}}
+ + ( \vec u^\prime \cdot \nabla^\prime ) \vec u^\prime =
+-\nabla^\prime p^\prime {\mu \over {\rho \kappa}} +
+\delta\rho \vec g^\prime {\mu \over {\rho \kappa}} +
+{\mu\over\rho} \left (\nabla^\prime \right )^2 \vec u^\prime {1\over {kappa}}.
+\label{eq:navier-stokes-ndpr}
+\end{equation}
+We can simplify Equation~\ref{eq:navier-stokes-ndpr} by collecting terms
+\begin{equation}
+{{\partial \vec u^\prime} \over {\partial t^\prime}} + ( \vec u^\prime \cdot
+\nabla^\prime ) \vec u^\prime = -\nabla^\prime p^\prime {\rm Pr}
++ \delta\rho \vec g^\prime {\rm Pr} +
+{\rm Pr} \left (\nabla^\prime \right )^2 \vec u^\prime.
+\label{eq:navier-stokes-nd2pr}
+\end{equation}
+where
+\begin{equation}
+Pr = {\mu\over{\rho\kappa}} = {\nu \over \kappa}
+\end{equation}
+is the Prandtl number. Prandtl was the director of the Kaiser Wihelm Institute
+for Fluid Research at the University of G\"ottingen. His student, Tiejens,
+wrote a classic text, ``Fundamentals of Hydro- and Aeromechanics'' (1934) which
+is based on Prandtl's lecture notes. It is an interesting and important
+attempt to bridge the gap between the theoretical fluid dynamics world, where
+approximations, even if not valid, were made to solve equations and the
+experimental world which was a collection of disintegrated and
+(seemingly) unrelated problems.
+
+Once again, the Prandtl number, like the Reynolds number, is dimensionless. Lets
+think about the physical meaning of the Prandtl number. One way you can think
+of it is the ratio of the diffusion of momentum, the dynamic viscosity,
+$\nu,$ to the diffusion of heat, $\kappa$, the thermal diffusivity. Notice that
+if I let the Prandlt number go to infinity we recover Stokes equation. The
+table below gives you the Prandtl number of some common fluids.
+
+\begin{center}
+\begin{tabular}{||c|c|c|c|c||} \hline
+Fluid & Density & Viscosity & Thermal Diff- & Prandtl number \\
+ & kg/m$^3$ & Pa s & usivity (m$^2$/s) & (dimensionless)\\ \hline
+Air & 1.3 & $10^{-5}$ & $2\times 10^{-5}$ & 0.72 \\
+Water & 1000.0 & $10^{-3}$ &$1.4\times 10^{-7}$& 8.10 \\
+Olive Oil & 916. & 0.1 &$9.2\times 10^{-8}$& 1,170 \\
+Glycerine & 1260. & 2.3 &$9.8\times 10^{-8}$& 18,880 \\
+Mantle & 3300.0 & $10^{21}$ & $10^{-6}$ & $10^{30}$ \\ \hline
+\end{tabular}
+\end{center}
+
+In practice, experimentalists usually consider a Prandtl number on the order of
+1,000 to be large enough to use the infinite Prandtl number approximation.
+
+\bigskip
+To convince yourselves that you understand non-dimensionalization, attempt the
+following problems.
+
+\begin{quote} {\bf Problem 3} Return to Equation~\ref{eq:navier-stokes1} and use
+the following relations between the dimensional and non-dimensional variables
+\begin{eqnarray}
+x & = & x^\prime D \\
+u & = & u^\prime {V} \\
+t & = & t^\prime {D\over V} \\
+p & = & P_o + p^\prime {\rho V^2} \\
+\rho g & = & \rho {V^2 \over D}
+\end{eqnarray}
+What is the dimensionless form of the equation?
+\end{quote}
+
+\begin{quote} {\bf Problem 4} You are asked to design a lab experiment to study
+flow in the Wabash River. Given that your lab space and budget will allow
+you to build a tank that is 2 meters deep and you will use water, what
+velocities will you need in the lab to reproduce the conditions on the Wabash
+(where the water flows at between 1-10 meters/second). Defend your answer.
+\end{quote}
+
+\vfill\eject
+
+\subsection{Solution Strategies}
+
+\subsubsection{One-Dimensional Channel Flow}
+
+We now turn to some simple solutions for low Reynolds number, or infinite
+Prandtl number, or viscous flow. For the first problem, we will restrict our
+attention to steady flow in one dimension. The figure blow illustrates the
+geometry of this flow.
+
+\begin{figure}[h]
+\begin{picture}(300,150)(0,0)
+\put(12,25){\framebox(275,100){}}
+\end{picture}
+\caption{The 1-D Channel}
+\end{figure}
+
+We will assume the fluid is a constant viscosity, Newtonian fluid. Thus we can
+use,
+\begin{equation}
+{\mu} \nabla^2 \vec u = \nabla p + \delta\rho \vec g.
+\end{equation}
+We will assume that $\mu$ and $\rho$ are constant and that $ \delta\rho = 0$,
+so there are no buoyancy forces driving this problem. Remember that we can
+write the equations for each component separately,
+\begin{equation}
+\mu \left ( {{\partial^2 u}\over{\partial x^2 }} + {{\partial^2
+u}\over{\partial z^2}} \right ) = {{\partial P}\over {\partial x}}.
+\end{equation}
+and
+\begin{equation}
+\mu \left ( {{\partial^2 v}\over{\partial x^2 }} + {{\partial^2
+v}\over{\partial z^2}} \right ) = {{\partial P}\over {\partial z}}.
+\end{equation}
+Because the flow is only in the $x$ direction, we can ignore the second
+equation ($v = 0$). Notice that this also tells us that there are no
+pressure gradients in the $z$ direction (${{\partial P}\over {\partial z}} =
+0$). In addition, we will assume ${{\partial u}\over {\partial x}} = 0$.
+However, we can not assume that ${{\partial P}\over {\partial x}} =0$, and
+could drive flow in the $u$ direction.
+
+If you wish, you could balance the forces on the box, or you can return to our
+equations above and get,
+\begin{equation}
+\label{eq:1d}
+\mu {{\partial^2 u}\over{\partial z^2}} = {{\partial P}\over {\partial x}}.
+\end{equation}
+
+If we assume ${{\partial P}\over {\partial x}}$ is independent of $z$ (which it
+is, ask yourself why) we can integrate equation~\ref{eq:1d} twice and get
+\begin{equation}
+u = {1\over{2\mu}} {{\partial P}\over {\partial x}} z^2 + c_1 z + c_2
+\end{equation}
+where $c_1$ and $c_2$ arise because we have used indefinite integrals.
+
+We can solve this simple problem for specific cases, by choosing boundary
+conditions which allow us to solve for $c_1$ and $c_2$.
+
+For our first example, lets consider ${{\partial P}\over {\partial x}} =0$ and
+a channel of height $z=D$. Lets choose $u(0) = 0$ and $u(D) = U_0$. From
+$u(0) = 0$ we immediately see that $c_2 = 0$. From $u(D) = U_0$, we find that
+$c_1 = {U_0 \over D}$. Thus the solution is
+\begin{equation}
+u(z) = z{U_0 \over D}.
+\end{equation}
+There are several things of interest here, first, the velocity is linear in $z$
+and second, the velocity is independent of the viscosity $\mu$. In this case,
+the boundary conditions completely determine the solution (independent of the
+material properties). This is called {\em Couette flow}.
+\begin{figure}[h]
+\begin{picture}(300,150)(0,0)
+\put(12,25){\framebox(275,100){}}
+\end{picture}
+\caption{Couette Flow}
+\end{figure}
+
+For our second example, lets consider the case where ${{\partial P}\over
+{\partial x}} \neq 0$. Lets eliminate the effect of non-zero velocity boundary
+conditions, which we looked at in the example above and consider, $u(0) =
+u(D) = 0$. From $u(0) = 0$, we see that $c_2$ = 0. From $ u(D) = 0$, we find
+\begin{equation}
+0 = {1\over{2\mu}} {{\partial P}\over {\partial x}} D^2 + c_1 D
+\end{equation}
+or
+\begin{equation}
+c_1 = {-1\over{2\mu}} {{\partial P}\over {\partial x}} D
+\end{equation}
+thus,
+\begin{equation}
+u(z) = {1\over{2\mu}} {{\partial P}\over {\partial x}} \left ( z^2 - Dz \right
+)
+\end{equation}
+
+Notice in this case that the solution is parabolic and by taking the
+derivative of $u(z)$ and setting it equal to zero, we find the maximum
+(minimum) is at $D/2$. This is called {\em Poiseuille Flow.}
+\begin{figure}[h]
+\begin{picture}(300,150)(0,0)
+\put(12,25){\framebox(275,100){}}
+\end{picture}
+\caption{Poiseuille Flow}
+\end{figure}
+
+A note here, because the equations of viscous flow are linear, we can
+superimpose two solutions. As long as we satisfy the boundary conditions
+correctly, they the superposition of two solutions is also a solution to the
+equation. Because we now have a solution with a non-zero boundary condition
+and a solution with a non-zero pressure gradient, we now have the tools to
+solve all the possible combinations of 1-D problems.
+
+We can, and people have, used these equations to gain intuition about viscous
+flow for Earth problems. A simple solution which we looked at last class was
+a simple model for plate motion and flow in the mantle. This allowed us to
+superimpose two solutions. The general solution, when we have both a velocity
+b.c. on one side of the channel and a pressure gradient is
+\begin{equation}
+u = {1\over{2\mu}} {{\partial P}\over {\partial x}}
+\left ( z^2 - zD \right ) - z{U_0 \over D} + U_0
+\end{equation}
+
+(The subsection on the asthenosphere counter-flow still needs to be typeset.
+Sorry).
+
+
+\subsubsection{Stream-Function Formulation}
+
+There are a number of strategies for solving viscous flow problems in more than
+one dimension, some of these can be done analytically (in special restricted
+cases), some can simplify the equations so that standard computer software or
+algorithms can be used. The first strategy we will look at is the
+stream-function technique. Lets return to the Stokes equation, with a constant
+viscosity.
+\begin{equation}
+\nabla^2 \vec u = \nabla p +{\bf Ra \,}T\hat k. \label{eq:const-stokes}
+\end{equation}
+This, along with Equation~\ref{eq:incompress-cont} can be solved in 2-D
+by using the stream-function. If we define a function $\psi$ such that
+\begin{eqnarray}
+u & = &-{{\partial \psi} \over {\partial y}} \\
+v & = & {{\partial \psi} \over {\partial x}}
+\end{eqnarray}
+where $u$ and $v$ are the $x$ and $y$ components of $\vec u$.
+Notice that Equation~\ref{eq:incompress-cont} is satisfied by our choice of
+$\psi $. Substituting our definitions of $u$ and $v$ into
+Equation~\ref{eq:const-stokes}, we get
+\begin{eqnarray}
+0 & = & {{\partial P} \over {\partial x}} +
+\mu \left ({{\partial^3 \psi} \over {\partial x^2 \partial y}} + {{\partial^3
+\psi} \over {\partial y^3}} \right ) \label{eq:x1} \\
+-{\bf Ra \,}T & = & -{{\partial P} \over {\partial y}} +
+\mu \left ({{\partial^3 \psi} \over {\partial x^3 }} + {{\partial^3\psi} \over
+{\partial x \partial y^2}} \right ) \label{eq:y1}
+\end{eqnarray}
+Taking the derivative with respect to $y$ of Equation~\ref{eq:x1}
+and the derivative with respect to $x$ of Equation~\ref{eq:y1} and adding them
+we get a biharmonic equation for the stream function $\psi , $
+
+\begin{equation}
+-{\bf Ra \,} {{\partial T}\over{\partial x}} =
+{{\partial^4 \psi} \over {\partial x^4}} +
+2 {{\partial^4 \psi} \over {\partial x^2\partial y^2}} +
+{{\partial^4 \psi} \over {\partial y^4}}. \label{eq:biharmonic}
+\end{equation}
+In in a more compact notation, this can be written as
+\begin{equation}
+-{\bf Ra \,} {{\partial T}\over{\partial x}} = \nabla^4 \psi
+\end{equation}
+The solution to the biharmonic equation is relatively straight-forward.
+
+Before we do that, lets examine the meaning of the stream function a bit more.
+We can write the stream function for the general one-dimensional channel flow
+problem that we did last time.
+
+The general solution, when we have a velocity b.c. on one side of the channel
+and a pressure gradient is
+\begin{equation}
+u = {1\over{2\mu}} {{\partial P}\over {\partial x}} \left ( z^2 - Zed \right ) -
+z{U_0 \over D} + U_0
+\end{equation}
+using the definition of the stream function,
+\begin{equation}
+u = -{{\partial \psi} \over {\partial y}}
+\end{equation}
+and integrating to find $\phi$, we get
+\begin{equation}
+\psi = {-1\over{2\mu}} {{\partial P}\over {\partial x}} \left (
+{{z^3}\over 3} - {{z^2 D}\over 2} \right ) + z^2{U_0 \over {2D}} - U_0
+\end{equation}
+
+It is {\bf not} helpful to plot $\psi$ for the Couette and Poiseuille Flows.
+
+Once again, to convince yourselves that you understand the 1-D problems, I
+suggest the following two problems:
+
+\begin{quote} {\bf Problem 5}
+Consider the steady, 1-D flow of a viscous fluid down an incline plane. Assume
+the flow occurs in a layer of thickness $h$, as shown in the figure below.
+Show that the velocity profile is given by:
+\begin{equation}
+u(y) = {{\rho g \sin \alpha}\over {2 \mu}} \left ( h^2 - y^2 \right )
+\end{equation}
+where $y$ is the coordinate measured perpendicular to the incline plane,
+$\alpha$ is the angle of inclination, and $g$ is the acceleration due to
+gravity.
+
+Hint: First show that
+\begin{equation}
+{{\partial \tau}\over {\partial y}} = - \rho g \sin \alpha
+\end{equation}
+then apply the no slip condition at $y=h$ and $\tau = 0$ at $y=0$. This is the
+{\em free-slip} condition.
+\end{quote}
+\begin{figure}[h]
+\begin{picture}(300,150)(0,0)
+\put(12,25){\framebox(275,100){}}
+\end{picture}
+\caption{Inclined Channel For Problem 5}
+\end{figure}
+
+\begin{quote} {\bf Problem 6} Suppose you have a moving plate sitting on top
+of a fluid that has two layers as pictured. Our solution for Couette flow
+only told us how to solve for a single viscosity fluid, but by matching the
+stresses just above and below the viscosity interface, you can solve this
+problem. Using the values in the picture, solve for the steady velocity
+profile.
+\end{quote}
+\begin{figure}[h]
+\begin{picture}(300,150)(0,0)
+\put(12,25){\framebox(275,100){}}
+\end{picture}
+\caption{Two Layer Channel For Problem 6}
+\end{figure}
+
+\subsection{Stream-Function Examples}
+
+\subsubsection{Post Glacial Rebound}
+
+The growth and melting of large ice sheets occurs on a time scale so fast
+compared to the mantle that loading and unloading the ice sheet can change flow
+in the mantle. During the last great ice age, Scandinavia and North America
+were covered with thick sheets of ice. This caused a considerable subsidence
+of the surface. When the ice sheets melting, the surface began to return to
+its normal elevation (rebound). The rate of rebound has been dated by changed
+in elevation of ancient shore lines.
+
+To study this problem, we will use the flow of a semi-infinite, viscous, half
+space. We will consider $y > 0$ subject to a periodic initial surface
+displacement given by
+\begin{equation}
+w_m = w_{m0} \cos {{2 \pi x}\over{\lambda}}
+\end{equation}
+where $\lambda$ is the wavelength of the initial displacement and $w_m \ll
+\lambda$. The displacement leads to a hydrostatic pressure gradient which
+drives flow, attempting to restore the surface to the undeformed state ($w =
+0$). We will use the stream-function formulation, solving the bi-harmonic
+equation, to determine the viscous flow. Because the initial displacement
+varies in $x$ like $\cos {{2 \pi x}\over{\lambda}}$ it is reasonable to assume
+that the stream function will also vary like a sin or cos in $x$. Recall that
+$\psi$ is an integral of velocity, so it is reasonable to expect that $\psi
+\propto \sin {{2 \pi x}\over{\lambda}}$. We will assume this form for the
+solution and assume that the depth variation can be separated from the
+horizontal variation, i.e.,
+\begin{equation}
+\psi = \sin {{2 \pi x}\over{\lambda}} Y(y)
+\end{equation}
+where $Y$ is the function we will solve for. If we substitute this form into
+the biharmonic equation~\ref{eq:biharmonic}, we get
+\begin{equation}
+\frac{d^4 Y}{dy^4} - 2 \left [ {{2 \pi}\over{\lambda}} \right ]^2
+\frac{d^2 Y}{dy^2} + \left [ {{2 \pi}\over{\lambda}} \right ]^4 Y = 0.
+\label{eq:forth-order-eq}
+\end{equation}
+Solutions of this equation will be of the form
+\begin{equation}
+Y(y) \propto \exp (my).
+\end{equation}
+If this form is substituted into the above equation for $Y(y)$, we find that
+$m$ is the solution of
+\begin{equation}
+m^4 - 2 \left [ {{2 \pi}\over{\lambda}} \right ]^2 m^2 + \left [ {{2
+\pi}\over{\lambda}} \right ]^4 = \left ( m^2 -\left [ {{2 \pi}\over{\lambda}}
+\right ]^2 \right )^2 = 0
+\end{equation}
+or $m = \pm {{2 \pi}\over{\lambda}}$. This provides for two of the four
+possible solutions for $Y$. By direct substitution, you can verify that
+\begin{eqnarray}
+Y(y) & \propto & y \exp \left [ {{2 \pi y}\over{\lambda}} \right ] \\
+Y(y) & \propto & y \exp \left [ -{{2 \pi y}\over{\lambda}} \right ]
+\end{eqnarray}
+are also solutions to~\ref{eq:forth-order-eq}. So the general solution can be
+written as
+\begin{equation}
+\psi = \sin \left (\frac{2 \pi x}{\lambda} \right ) \left [ A e^{\frac{-2 \pi
+y}{\lambda}} + By e^{\frac{-2 \pi y}{\lambda}} + C e^{\frac{2 \pi y}{\lambda}}
++ Dy e^{\frac{2 \pi y}{\lambda}} \right ]
+\end{equation}
+with the four constants $A, B, C, D$ to be determined by the boundary
+conditions. We first use the fact that we want finite velocities as the depth
+of the half-space goes to infinity. Therefore, $C$ and $D$ must be zero. So
+the formula for the stream function reduces to
+\begin{equation}
+\psi = \sin \left (\frac{2 \pi x}{\lambda} \right )
+e^{\frac{-2 \pi y}{\lambda}} \left [ A + By \right ]
+\end{equation}
+We can now differentiate $\psi$ to get the velocity components
+\begin{eqnarray}
+u & = & \sin \left (\frac{2 \pi x}{\lambda} \right )
+e^{\frac{-2\pi y}{\lambda}}\left [\frac{2\pi}{\lambda}(A + By)-B\right ] \\
+v & = & \cos \left (\frac{2 \pi x}{\lambda} \right ) \frac{2 \pi}{\lambda}
+e^{\frac{-2 \pi y}{\lambda}} \left [ A + By \right ]
+\end{eqnarray}
+Because of the extreme different in properties between the mantle and the
+atmosphere (here I am assuming that the lithosphere/crust is the same as the
+mantle, which for the wavelengths we are talking about is not too bad an
+assumption) we force the horizontal velocity of the flow to be zero at $y =
+w$. Now we use the fact that $w \ll \lambda$ so that we can apply this
+condition to $y=0$. By setting $u=0$ at $y=0$ we find that $B = \frac{2\pi
+A}{\lambda}$ which means
+\begin{eqnarray}
+\psi & = & A \sin \left (\frac{2 \pi x}{\lambda} \right ) e^{\frac{-2 \pi
+y}{\lambda}} \left [ 1 + \frac{2 \pi y}{\lambda} \right ] \\
+u & = & A \left ( \frac{2\pi}{\lambda} \right )^2
+\sin \left (\frac{2 \pi x}{\lambda} \right ) y e^{\frac{-2\pi y}{\lambda}} \\
+v & = & A \left ( \frac{2 \pi}{\lambda} \right )
+\cos \left (\frac{2 \pi x}{\lambda} \right ) e^{\frac{-2 \pi y}{\lambda}}
+\left [ 1 + \frac{2 \pi y}{\lambda} \right ]
+\end{eqnarray}
+In order to find the constant $A$, we need to equate the hydrostatic pressure
+head associated with the topography $w$ to the normal stress at the upper
+boundary. They hydrostatic pressure is $-\rho g w$ and the normal stress is
+$p - 2\mu \frac{\partial v}{\partial y}$. Once again because the
+displacements are small, this can be done at $y=0$. First, we must calculate
+the pressure at $y=0$, which we can do by putting our expressions for $u$ and
+$v$ back into the viscous flow equation (not the biharmonic). Doing this, we
+obtain
+\begin{equation}
+\frac{\partial p}{\partial x} = -2 \mu A \left ( \frac{2 \pi}{\lambda} \right
+)^3 \sin \frac{2 \pi x}{\lambda}
+\end{equation}
+Integrating this we get
+\begin{equation}
+p = 2 \mu A \left ( \frac{2 \pi}{\lambda} \right
+)^2 \cos \frac{2 \pi x}{\lambda}
+\end{equation}
+on $y=0$. We also need $\frac{\partial v}{\partial y}$ at $y=0$ which
+fortunately is $0$. Thus
+\begin{eqnarray}
+-\rho g w & = & p - 2\mu \frac{\partial v}{\partial y}\\
+w_{y=0} & = & \frac{-2 \mu A}{\rho g} \left ( \frac{2 \pi}{\lambda} \right )^2
+\cos\frac{2 \pi x}{\lambda}
+\end{eqnarray}
+
+The key step is to notice that the time derivative of displacement $w$ is the
+vertical velocity (at $y=w$). But once again we use the fact that $w \ll
+\lambda$ so that we can assume this at $y=0$. At $y-0$ we have
+\begin{equation}
+v_{y=0} = A \frac{2 \pi}{\lambda}
+\cos \left (\frac{2 \pi x}{\lambda} \right )
+\end{equation}
+so
+\begin{equation}
+\frac{d w}{dt}_{y=0} = A \frac{2 \pi}{\lambda}
+\cos \left (\frac{2 \pi x}{\lambda} \right )
+\end{equation}
+or
+\begin{equation}
+\frac{d w}{dt}_{y=0} = -w \frac{\lambda g \rho}{4 \pi \mu}
+\end{equation}
+we can integrate this to get
+\begin{equation}
+w(t) = w_{m0} \exp \left ( \frac{\lambda g \rho t}{4 \pi \mu} \right )
+\end{equation}
+thus the surface rebounds exponentially with time and fluid flows from regions
+of elevated topography to regions of low topography. We can think of the
+grouping of constants $ \frac{4 \pi \mu}{\lambda g \rho} $ as a characteristic
+time (i.e., the time it takes the topography to decay by 1/e). Notice that
+this grouping has units of time. If we know the density (we do approximately)
+and we know $g$ and the initial wavelength, $\lambda$, which we know pretty
+well, then we can estimate the viscosity.
+
+This exercise was first discussed by Norman Haskell, in 1935. Haskell came up
+with a value (in todays units) of $10^{21}$ Pa s which is sometimes called the
+Haskell value.
+
+\medskip
+\noindent
+Haskell, N.A., The motion of a viscous fluid under a surface load, I,
+{\em Physics,} 6, 265-269, 1935.
+
+\medskip
+\noindent
+Haskell, N.A., The motion of a viscous fluid under a surface load, II,
+{\em Physics,} 7, 56-61, 1936.
+
+It is important to note that this is an average value, or effective value if
+the Earth were a homogeneous fluid (which it is not). Subsequent work has
+focused on layered, then radial models (the spherical nature is not a
+significant factor for all the additional work it contributes). Layered
+half-spaces are a fairly easy step from the discussion here. It is also
+straight-forward to consider an elastic lithosphere over a viscous half-space.
+A detailed discussion of these effects can be found in
+
+\medskip
+\noindent
+Cathles, L.M., {\em The Viscosity of the Earth's Mantle,} 386 pp., Princeton
+University Press, Princeton, NJ, 1975.
+
+\medskip
+
+Current research focuses on laterally varying material properties. The
+methods that we will go on to discuss are needed to solve that level of problem.
+
+\subsubsection{Corner Flow and the Angle of Subduction}
+
+One of the most curious problems in mantle convection is that slabs do not
+descend vertically. All models of convection of a viscous fluid produce
+symmetric, two-sided downwelling unless something is added to break the
+symmetry, like a dipping fault or some difference in the properties of the two
+thermal boundary layers. Some data
+
+\begin{center}
+\begin{tabular}{||l|r||} \hline
+Arc & Dip Angle\\ \hline
+Central Chile & 5$^\circ$ \\
+Northern Chile & 30$^\circ$ \\
+Southern Chile & 30$^\circ$ \\
+Honshu & 30$^\circ$ \\
+Izu-Bonin & 60$^\circ$ \\
+Java & 70$^\circ$ \\
+New Hebrides & 70$^\circ$ \\
+Ryukyu & 45$^\circ$ \\
+West Indies & 50$^\circ$ \\ \hline
+\end{tabular}
+\end{center}
+
+One explanation for why the slab descends at an angle other than $90^\circ$ is
+that pressure forces due to the induced flow pull the slab up and
+counter-balance the buoyancy forces. We can investigate this with a stream
+function solution.
+
+As unlikely as it may seem, it is possible to write $\psi$ in the form
+\begin{equation}
+\psi ~=~ A\,x + B\,y + \left( {\rm C}\,x + {\rm D}\,y \right) \,
+ \arctan ({y\over x})
+\end{equation}
+where $A, B, C,$ and $D$ are constants to be determined by the boundary
+conditions. Because of the geometry of the descending slab, there are two
+distinct stream functions, i.e., two sets of constants to be solved for in this
+problem: one above and one below the slab.
+It is helpful to recall that
+\begin{eqnarray}
+\frac{\partial}{\partial y} \arctan \frac{y}{x} & = & \frac{x}{x^2+y^2} \\
+\frac{\partial}{\partial x} \arctan \frac{y}{x} & = & \frac{-y}{x^2+y^2}
+\end{eqnarray}
+Differentiating the stream function to find the velocities, we find that
+\begin{eqnarray}
+u(x,y) & = & -\frac{\partial\psi}{\partial y} ~=~ -B - {{{\rm C}\,x + {\rm
+D}\,y}\over{x\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }} -
+ {\rm D}\,\arctan ({y\over x}) \\
+v(x,y) & = & \frac{\partial\psi}{\partial x} ~=~ A - {{y\,\left( {\rm C}\,x +
+{\rm D}\,y \right) }\over{{x^2}\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }} +
+ {\rm C}\,\arctan ({y\over x}).
+\end{eqnarray}
+
+First, notice that this will satisfy continuity, because
+\begin{equation}
+\frac{\partial u}{\partial x} ~=~ {{-2\,{y^2}\,\left( {\rm C}\,x + {\rm D}\,y
+\right) }\over
+ {{x^4}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 2}}} - {{{\rm C}}\over
+ {x\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }} +
+ {{{\rm D}\,y}\over
+ {{x^2}\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }} +
+ {{{\rm C}\,x + {\rm D}\,y}\over
+ {{x^2}\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }} \label{eq:ux}
+\end{equation}
+and
+\begin{equation}
+\frac{\partial v}{\partial y} ~=~{{2\,{y^2}\,\left( {\rm C}\,x + {\rm D}\,y
+\right) }\over
+ {{x^4}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 2}}} + {{{\rm C}}\over
+ {x\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }} -
+ {{{\rm D}\,y}\over
+ {{x^2}\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }} -
+ {{{\rm C}\,x + {\rm D}\,y}\over
+ {{x^2}\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }} \label{eq:vy}
+\end{equation}
+By inspection ~\ref{eq:ux} plus \ref{eq:vy} sum to zero.
+
+
+As in the last problem, the pressure can be found by substituting the solution
+for $u$ back into the equation of motion and solving for $\frac{\partial
+p}{\partial x}$ and integrating in $x$. To see this notice that
+\begin{displaymath}
+\frac{\partial^2 u}{\partial x^2} ~=~ {{-8\,{y^4}\,\left( {\rm C}\,x + {\rm
+D}\,y \right) }\over
+ {{x^7}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 3}}} - {{4\,{\rm C}\,{y^2}}\over
+ {{x^4}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 2}}} + {{2\,{\rm D}\,{y^3}}\over
+ {{x^5}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 2}}} +
+\end{displaymath}
+\begin{equation}
+{{10\,{y^2}\,
+ \left( {\rm C}\,x + {\rm D}\,y \right) }\over
+ {{x^5}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 2}}} + {{2\,{\rm C}}\over
+ {{x^2}\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }} -
+ {{2\,{\rm D}\,y}\over
+ {{x^3}\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }} -
+ {{2\,\left( {\rm C}\,x + {\rm D}\,y \right) }\over
+ {{x^3}\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }}
+\end{equation}
+and
+\begin{equation}
+\frac{\partial^2 u}{\partial y^2} ~=~{{-8\,{y^2}\,\left( {\rm C}\,x + {\rm D}\,y
+\right) }\over
+ {{x^5}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 3}}} + {{6\,{\rm D}\,y}\over
+ {{x^3}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 2}}} + {{2\,\left( {\rm C}\,x + {\rm D}\,y \right)}
+ \over {{x^3}\,{{\left( 1 + {{{y^2}} \over {{x^2}}} \right
+ ) }^2}}}
+\end{equation}
+thus,
+\begin{eqnarray}
+\frac{\partial p}{\partial x} & = &\mu \left ( \frac{\partial^2 u}{\partial x^2}
++ \frac{\partial^2 u}{\partial y^2} \right ) = \nonumber \\
+&& \mu \,\left( {{-8\,{y^2}\,
+ \left( {\rm C}\,x + {\rm D}\,y \right) }\over
+ {{x^5}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 3}}} - {{8\,{y^4}\,
+ \left( {\rm C}\,x + {\rm D}\,y \right) }\over
+ {{x^7}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 3}}} + {{6\,{\rm D}\,y}\over
+ {{x^3}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 2}}} - \right. \nonumber \\
+ && {{4\,{\rm C}\,{y^2}}\over
+ {{x^4}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 2}}} + {{2\,{\rm D}\,{y^3}}\over
+ {{x^5}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 2}}} + {{2\,\left( {\rm C}\,x + {\rm D}\,y \right )
+ }\over
+ {{x^3}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right ) }^
+ 2}}} + {{10\,{y^2}\,
+ \left( {\rm C}\,x + {\rm D}\,y \right) }\over
+ {{x^5}\,{{\left( 1 + {{{y^2}}\over {{x^2}}} \right) }^
+ 2}}} + \nonumber \\
+ && \left. {{2\,{\rm C}}\over
+ {{x^2}\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }} -
+ {{2\,{\rm D}\,y}\over
+ {{x^3}\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }} -
+ {{2\,\left( {\rm C}\,x + {\rm D}\,y \right) }\over
+ {{x^3}\,\left( 1 + {{{y^2}}\over {{x^2}}} \right) }} \right )
+\label{eq:harry-mess}
+\end{eqnarray}
+and finally integrating~\ref{eq:harry-mess} with respect to $x$ gives
+\begin{equation}
+p(x,y) ~=~ {{-2\,\mu \,\left( {\rm C}\,x + {\rm D}\,y \right) }\over
+ {{x^2} + {y^2}}}
+\end{equation}
+Now we are ready to apply boundary conditions to the problem.
+
+The general problem leads to quite complicated integrals. We can look at some
+specific cases. A dip of $\pi/4$ is a good choice. We will first look at the
+solution in the arc corner. We have the following boundary conditions
+\begin{eqnarray}
+u=v & = & 0 ~~~~~~~~~~~~~~~~~~~~~on~y=0,~x>0~or~\arctan \frac{y}{x}~=~0 \\
+u=v & = & {U\sqrt{2}\over 2}~~~~~~~~~~~~~~~on~y=x~or~\arctan
+\frac{y}{x}~=~\frac{\pi}{4}
+\end{eqnarray}
+Substituting $y=0$ into $u(x,y) = 0$, we find
+\begin{equation}
+-B -C = 0. \label{eq:bc1}
+\end{equation}
+Substituting $y=0$ into $v(x,y) = 0$, we find
+\begin{equation}
+A=0. \label{eq:bc2}
+\end{equation}
+Substituting $y=x (\arctan\frac{y}{x} = \pi/4)$, into $u(x,y) = U\sqrt{2}/2$, we
+find
+\begin{equation}
+-B - \frac{C}{2} -\frac{D}{2} -\frac{D\pi}{4} = U\sqrt{2}/2. \label{eq:bc3}
+\end{equation}
+Substituting $y=x (\arctan\frac{y}{x} = \pi/4)$, into $v(x,y) = U\sqrt{2}/2$, we
+find
+\begin{equation}
+A - \frac{C}{2} -\frac{D}{2} + \frac{C\pi}{4} = U\sqrt{2}/2.
+\label{eq:bc4}
+\end{equation}
+Simultaneously solving~\ref{eq:bc1}, \ref{eq:bc2}, \ref{eq:bc3}, and
+\ref{eq:bc4}, leads to the following expressions for
+$C$ and $D$,
+\begin{eqnarray}
+C & = & \frac{-\pi U\sqrt{2}}{2(2-\pi^2/4)} \\
+D & = & \frac{-U\sqrt{2}(2-\pi/2)}{(2-\pi^2/4)}
+\end{eqnarray}
+thus the pressure in the arc corner is
+\begin{equation}
+p = \frac{\mu U \sqrt{2} \left [ \pi x + (4-\pi)y \right]}{(2-\pi^2/4)(x^2+y^2)}
+\end{equation}
+we can evaluate this expression using the fact that $x=y=r\sqrt{2}/2$ where $r$
+is the length along the slab. Thus we find that
+\begin{equation}
+p = \frac{4\mu U}{(2-\pi^2/4)r} = \frac{-8\mu U}{r}
+\end{equation}
+The negative pressure means the back-arc region flow tends to lift the slab
+against the force of gravity. Notice that the pressure varies as $1/r$ along
+the slab, with a singularity at the trench. The total lifting torque on the
+slab will be $r\times P$ which is constant.
+
+For the ocean side of the arc, the boundary conditions are
+\begin{eqnarray}
+u=U, v & = & 0 ~~~~~~~~~~~~~~~~~~~on~y=0,~x<0~or~\arctan \frac{y}{x}~=~\pi \\
+u=v & = & {U\sqrt{2}\over 2}~~~~~~~~~~~~~on~y=x~or~\arctan
+\frac{y}{x}~=~\frac{\pi}{4}
+\end{eqnarray}
+Substituting $y=0, (\arctan \frac{y}{x}= \pi )$ into $u(x,y) = U $ we get
+\begin{equation}
+-B - C - D \pi = U. \label{eq:bc5}
+\end{equation}
+Substituting $y=0, (\arctan \frac{y}{x}= \pi )$ into $v(x,y) = 0 $ we get
+\begin{equation}
+A + C \pi = 0. \label{eq:bc6}
+\end{equation}
+Substituting $y=x, (\arctan \frac{y}{x}= \pi/4 )$ into $u(x,y) = U\sqrt{2}/2 $ we
+get
+\begin{equation}
+-B - \frac{C}{2} -\frac{D}{2} -\frac{D\pi}{4} = U\sqrt{2}/2. \label{eq:bc7}
+\end{equation}
+Substituting $y=x, (\arctan \frac{y}{x}= \pi/4 )$ into $v(x,y) = U\sqrt{2}/2 $
+we get
+\begin{equation}
+A - \frac{C}{2} -\frac{D}{2} + \frac{C\pi}{4} = U\sqrt{2}/2. \label{eq:bc8}
+\end{equation}
+Solving \ref{eq:bc5}, \ref{eq:bc6}, \ref{eq:bc7}, and \ref{eq:bc8} for $C$ and
+$D$, we get
+\begin{eqnarray}
+C & = & \frac{U}{(9\pi^2/4-2)} \left [ 2 - \frac{\sqrt{2}}{(1-3\pi/2)} \left (
+\frac{3\pi}{2} + \frac{9\pi^2}{4} \right ) \right ] \\
+D & = & \frac{U}{(9\pi^2/4-2)} \left [\sqrt{2} \left (2+\frac{3\pi}{2} \right )
+-2 \left ( 1+\frac{3\pi}{2} \right ) \right ]
+\end{eqnarray}
+The pressure on the bottom of the slab is given by
+\begin{equation}
+p= \frac{\mu U}{r} \left ( \frac{3\pi\sqrt{2} -4}{9\pi^2/4-2} \right ) =
+\frac{0.5 \mu U}{r}.
+\end{equation}
+Thinking about the geometry and the signs of the pressures, this flow also
+exerts a lifting torque on the slab, it is about a factor of 20 smaller than
+the suction torque. However, because constant viscosity convection never
+produces this kind of solution, we are lead to suspect that the downward force
+of gravity far outweighs the pressure torques. It would be an interesting
+exercise to see if the problem with both sides coming together at the same
+velocity produced torques that balanced at some angle...think about this.
+
+\end{document}
\ No newline at end of file
Added: mc/2D/ConMan/trunk/doc/viscous-flow2b.tex
===================================================================
--- mc/2D/ConMan/trunk/doc/viscous-flow2b.tex (rev 0)
+++ mc/2D/ConMan/trunk/doc/viscous-flow2b.tex 2008-06-30 17:59:49 UTC (rev 12357)
@@ -0,0 +1,1427 @@
+\documentclass{article}
+\usepackage{times}
+\begin{document}
+
+
+\title{EAS 591D Finite Element Methods for Viscous Flow}
+\author{Scott D. King}
+\date{Spring 2006}
+\maketitle
+
+\pagestyle{myheadings}
+\markright{Finite Element Methods for Viscous Flow}
+\setcounter{section}{1}
+\setcounter{page}{40}
+\section{Finite Element Method}
+
+Now we begin our introduction to the finite element method. To start with, I
+will consider the Possion Equation, because it avoids a complication that we
+will come to in a few weeks, while illustrating all the important points about
+second order differential equations. Recall Possion's equation
+\begin{equation}
+\nabla^2 u = f.
+\end{equation}
+
+I am going to follow a similar form to Hughes (The Finite Element Method)
+Chapter~1 verses 1-4. Lets think today in 1-D. There are two forms we can
+write the equation, the strong and the weak form. You are used to seeing the
+strong form, but not used to seeing the weak form. The finite element method,
+is cast in the weak form. In elasticity, for example, the weak form comes
+from a variational principal, such as the principal of virtual displacements.
+For viscous flow, there is also a variational form, but we will not go into
+that.
+
+\subsection{The Strong Form}
+
+The strong form can be formally stated: given $f: \bar \Omega \rightarrow \Re $
+and constants $g$ and $h$ find $u: \bar \Omega \rightarrow \Re$
+\begin{eqnarray}
+u,_{xx} ~+~ f & = & 0 \\
+u(1) & = & g \\
+u,_x (0) + h & = & 0
+\end{eqnarray}
+This choice of initial conditions allows us to examine both kinds of boundary
+conditions. The solution is trivial, but that does not matter. For
+completeness, it is
+\begin{equation}
+u(x) = g + (1-x) h + \int_x^1 \left ( \int_0^y f(z) dz \right ) dy
+\end{equation}
+
+
+\subsection{The Weak Form}
+
+To state the weak form formally, we define a set of trial solutions
+\begin{equation}
+\delta ~=~ \{ u~|~u ~\in H^1 , ~~u(1)=g \}
+\end{equation}
+and define a set of weighting functions
+\begin{equation}
+\varphi ~=~ \{ w~|~w ~\in H^1 , ~~w(1)=0 \}
+\end{equation}
+Given $f, g,$ and $h$ as before. Find $u \in \delta$ such that for all
+$w \in \varphi$
+\begin{equation}
+\int_0 ^1 w,_x u,_x dx ~=~ \int_0 ^1 w f dx + w(0) h
+\end{equation}
+In mechanics, the $w$'s are the virtual displacements, and the weak form is the
+variational equation.
+
+\subsection{Equivalence of Stong and Weak forms}
+
+We need to define a concept of a square integrable function, sometimes writen as
+$H^1$-functions. A function is square integrable on
+$[0,1]$ if
+\begin{equation}
+\int_0^1 (u,_x)^2 dx < \infty.
+\end{equation}
+
+We can show that the strong and weak forms of the equation are equivalent.
+Since both possess unique solutions (I will leave that for pure mathematicians
+to prove) the strong and weak solutions are one in the same.
+
+We begin with $u$ which we will assume is a solution to the strong form. Thus
+we can write
+\begin{equation}
+0 = - \int_0 ^1 w \{ u,_{xx} + f \} dx
+\end{equation}
+for $w \in \varphi$. Using integration by parts:
+\begin{equation}
+0 ~=~ \int_0 ^1 w,_x u,_x dx - \int_0 ^1 w\, f dx - w\, u,_x |^{^1}_{_0}
+\end{equation}
+so
+\begin{equation}
+\int_0 ^1 w,_x u,_x dx ~=~ \int_0 ^1 w\, f\, dx + w(0) h.
+\end{equation}
+This is the statement of the weak form of the problem. Furthermore, since $u$
+is a solution of the strong form $u(1) = g$ and is in $\delta$.
+
+We can work the problem backwards, starting with the weak form, using
+integration by parts and the definition of square integrability to prove the
+weak form is equivalent to the strong form. We begin with $v$ a solution to
+the weak form
+\begin{equation}
+\int_0 ^1 w,_x v,_x dx ~=~ \int_0 ^1 w\, f\, dx + w(0) h.
+\end{equation}
+We will use integration by parts and $w(1) = 0$ to get
+\begin{equation}
+0 = -\int_0 ^1 w \{ v,_{xx} + f \} dx + w(0) \{ v,_x (0) + h \}.
+\label{eq:intbyparts}
+\end{equation}
+To prove $v$ is a solution to the strong form, $ v,_{xx} + f = 0$ on
+$\Omega$ and $v,_x (0) + h = 0$. The first can be accomplished by choosing
+\begin{equation}
+w = \phi \left ( v,_{xx} + f \right )
+\end{equation}
+where $\phi$ is smooth, greater than zero inside the domain and zero at the
+boundaries. A good choice might be $\phi = x \left ( 1 -x \right )$.
+Subsitiuting our choice of $w$ into~\ref{eq:intbyparts} we get
+\begin{equation}
+0 = \int _0 ^ 1 \phi \left ( v,_{xx} + f \right )^2 dx + 0
+\end{equation}
+Because the squared term must be greater than or equal to zero, and $\phi >
+0$ in the domain, the only possible solution is that $v,_{xx} + f = 0$ within
+the domain $\Omega$. Using that fact, we can substitute back into
+Equation~\ref{eq:intbyparts} to get $v,_x (0) + h = 0$.
+
+\subsection{Galerkin's Approximation}
+
+Now we have a start on the finite element method. I want to continue to follow
+Hughes; however his notation becomes quite difficult to keep up with. Now,
+lets begin to think about putting a solution on the computer. Because we will
+have a finite approximation, related to how fine we space our grid, our
+solution will only approximate the real solution. Following Hughes' notation,
+the solution on the grid will be denoted as $u^h$ where $h$ is some measure of
+the spacing at the grid. Then,
+\begin{equation}
+\int_0 ^1 w^h,_x u^h,_x dx ~=~ \int_0 ^1 w^h\, f^h\, dx + w^h(0) h.
+\label{eq:weak}
+\end{equation}
+approximates our exact solution $u$.
+
+On a computer, we don't have a continuous solution. We have a solution at
+descrete points. We need to approximate the solution between the points (in
+order to integrate over the function). We will do this with {\bf shape
+functions} as they are usually called in the finite element language. Hughes
+uses $N_A ~~~A=1, 2, \cdots , n$ to denote the shape functions. You can also
+think of these as basis functions or interpolation functions. We require $ N_A
+(1) = 0, A=1,2, \cdots, n$. In order to specify our boundary condition, we
+need another shape function which has the property
+\begin{equation}
+N_{n+1} (1) = 1.
+\end{equation}
+Then, $g^h$ is given by,
+\begin{equation}
+g^h = g N_{n+1}
+\end{equation}
+and thus,
+\begin{equation}
+g^h (1) = g.
+\end{equation}
+With these definitions, we can write our solution $u^h$ as
+\begin{equation}
+u^h = \sum_{A=1}^{n} \, d_A \, N_A + g N_{n+1}
+\end{equation}
+where the $d_A$'s are unknown constants to be solved for.
+
+At this point, I don't want you to become discouraged, because we will make the
+shape functions more concrete. I want you to see how general this is, because
+in principle there is a great deal of flexibility in how we choose the shape
+functions.
+
+We have not said anything more about this function $w^h$ and how we are going
+to choose it. If our shape functions form a basis set for the grid, then we
+can represent {\bf any} function as a sum of the basis functions times some
+arbitrary coefficients $c_i$,
+\begin{equation}
+w^h = \sum_{A=1}^{n} \, c_A \, N_A = c_1 \, N_1 ~+~ c_2 \, N_2 ~+~ \cdots ~+~
+c_n \, N_n \label{eq:whapprox}
+\end{equation}
+If you don't remember this part of your mathematics background think of Fourier
+series. Any function one-dimensional function can be represented as an
+infinite series of sines and cosines times some unique set of coefficients.
+The shape functions form a similar kind of basis set.
+
+Notice that because we required that $ N_A (1) = 0, A=1,2, \cdots, n$,
+Equation~\ref{eq:whapprox} satisfies the requirement that $w^h (1) = 0$, as
+necessary.
+
+Using our definitions of the $w^h$'s and our approximation for $u^h$, we can
+get the messy expression for Equation~\ref{eq:weak}
+\begin{displaymath}
+\int_0 ^1 \left ( \frac{\partial}{\partial x}
+\left ( \sum_{A=1}^{n} \, c_A \, N_A \right )
+\frac{\partial}{\partial x}
+\left ( \sum_{B=1}^{n} \, d_B \, N_B + g N_{n+1} \right ) \right ) dx ~=~
+\end{displaymath}
+\begin{equation}
+\int_0 ^1 \sum_{A=1}^{n} \, c_A \, N_A \, f^h\, dx +
+\sum_{A=1}^{n} \, c_A N_A (0) h.
+\end{equation}
+By rearranging, we can write
+\begin{equation}
+\sum_{A=1}^{n} G_A c_A = 0
+\end{equation}
+where
+\begin{displaymath}
+G_A ~=~ \int_0 ^1 \left ( \frac{\partial N_A}{\partial x} \right )
+\left ( \sum_{B=1}^{n} \, d_B \, \frac{\partial N_B}{\partial x} \right ) dx
+\end{displaymath}
+\begin{equation}
+- \int_0 ^1 N_A \, f^h\, dx - N_A (0) h +
+\int_0 ^1 \frac{\partial N_A}{\partial x} \frac{\partial N_{n+1}}{\partial x}
+\, g\, dx
+\end{equation}
+Now I use the fact that the shape functions are basis functions, so
+$N_A \times N_B$ is zero except when $A = B$. We could equally well use the
+fact that the $c_A$'s are arbitrary. Both of these force us to conclude that
+each $G_A$ must be identically zero and we get
+\clearpage
+\begin{displaymath}
+\sum_{B=1}^{n} \, \left ( \int_0 ^1 \frac{\partial N_A}{\partial x}
+\frac{\partial N_B}{\partial x} \, dx \right ) \, d_B =
+\end{displaymath}
+\begin{equation}
+\int_0 ^1 N_A \, f^h\, dx + N_A (0) h -
+g\, \int_0 ^1 \frac{\partial N_A}{\partial x}
+\frac{\partial N_{n+1}}{\partial x} \, dx \label{eq:feform}
+\end{equation}
+Everything in Equation~\ref{eq:feform} is known except the $d_B$'s. This
+constitutes a system of $n$ equations and $n$ unknowns. We can think of the
+left hand side as a matrix, $K_{AB}$ whose entries are
+\begin{equation}
+\int_0 ^1 \frac{\partial N_A}{\partial x} \frac{\partial N_B}{\partial x} \, dx
+\label{eq:Kelements}
+\end{equation}
+We can write
+\begin{equation}
+\sum_{B=1}^{n} \, K_{AB} d_B = F_A, ~~~~~ A = 1, 2, \cdots, n
+\end{equation}
+or as a matrix equation
+\begin{equation}
+\left [ K \right ] \, \{d\} = \{f\}
+\end{equation}
+where
+\begin{equation}
+[k] = \left[
+\begin{array}{cccc}
+K_{11} & K_{12} & \cdots & K_{1n} \\
+K_{21} & K_{22} & \cdots & K_{2n} \\
+\vdots & \vdots & & \vdots \\
+K_{n1} & K_{n2} & \cdots & K_{nn}
+\end{array} \right ] \label{Kmatrix}
+\end{equation}
+By tradition, $\left [ K \right ]$ is the stiffness matrix, $\{f\}$ is the force
+vector, and $\{d\}$ is the displacement vector. When the problem under
+consideration pertains to a mechanical system, this makes the most sense, but
+even in heat conduction problems, or fluid flow problems, the terminology
+is still (often) retained.
+
+\subsection{Some Simple Problems}
+
+Lets look at a problem with one degree of freedom ($n\, =\, 1$). It will
+illustrate some important concepts. If $n\, =\, 1$, then $w^h = c_1 \, N_1$
+and $u^h \, = \, d_1 N_1 ~+~ g N_2$. The only unknown is $d_1$. The shape
+functions must satisfy $N_1(1) ~=~ 0$ and $N_2(1) ~=~ 1$. Lets choose
+\begin{eqnarray}
+N_1(x) ~=~ 1 \, - \, x \label{eq:n1} \\
+N_2(x) ~=~ x \label{eq:n2}
+\end{eqnarray}
+clearly $N_1$ and $N_2$ satisfy our required conditions. Because we have only
+one degree of freedom, the matrix collapses as follows:
+\begin{eqnarray}
+[K] ~=~ K_{11} \\
+\{f\} ~=~ F_1 \\
+\{d\} ~=~ d_1
+\end{eqnarray}
+Using Equation~\ref{eq:Kelements}, and substituting Equations~\ref{eq:n1}
+and~\ref{eq:n2} we find
+\begin{equation}
+K_{11} ~=~ \int _0 ^1 N_{1,x} N_{1,x} \, dx ~=~ \int _0 ^1 \{ -1\} \{-1\} \, dx
+~=~ 1
+\end{equation}
+and
+\begin{eqnarray}
+F_1 &~=~& \int _0 ^1 N_1(x) \, f(x) \, dx ~+~ N_1(0) \, h ~-~ g \, \int _0 ^1
+N_{1,x} \, N_{2,x} \, dx \nonumber \\
+ &~=~& \int _0 ^1 \left (1\, -\, x \right ) f(x) \, dx ~+~ h ~-~
+g \, \int _0 ^1 \{ -1 \} \{ 1 \} \, dx \nonumber \\
+ &~=~& \int _0 ^1 \left (1\, -\, x \right ) f(x) \, dx ~+~ h ~+~ g
+\end{eqnarray}
+so
+\begin{equation}
+d_1 ~=~ K_{11}^{-1} F_1 ~=~ F_1
+\end{equation}
+So if we substitute $d_1, N_1,$ and, $ N_2$ back into our expression for $u^h$,
+we get
+\begin{eqnarray}
+u^h \, = \, d_1 N_1 ~+~ g N_2 \\
+u^h (x) \, = \, \left ( \int _0 ^1 \left (1\, -\, y \right ) f(y) \, dy ~+~ h ~+~
+g \right ) \left ( 1 ~-~ x \right ) ~+~ g \, x \\
+\end{eqnarray}
+which is the finite element approximate solution for $u(x)$.
+
+It is helpful to compare this to the exact solution from last class,
+\begin{equation}
+u(x) ~=~ g ~+~ (1 - x) \, h ~+~ \int_x^1 \left ( \int_0^y f(z) dz \right ) dy
+\end{equation}
+
+Lets look as several specific cases:
+
+\noindent
+{\bf Case 1:} $f(x) = 0 $. In this case
+\begin{eqnarray}
+u(x) &=& g ~+~ (1 - x) \, h \\
+u^h (x) \, &=& (h + g) \, ( 1 ~-~ x ) ~+~ g \, x \nonumber \\
+ &=& g ~+~ (1 - x) \, h
+\end{eqnarray}
+In this case, $u(x) = u^h(x)$.
+
+\noindent {\bf Case 2:} $f(x) = p ~$ where $p$ is a constant. In this case
+\begin{eqnarray}
+u(x) &=& g ~+~ (1 - x) \, h ~+~ \frac{p(1-x^2)}{2}\\
+u^h (x) \, &=& (h + g) \, ( 1 ~-~ x ) ~+~ g \, x + \frac{p(1-x)}{2}\nonumber \\
+ &=& g ~+~ (1 - x) \, h + \frac{p(1-x)}{2}
+\end{eqnarray}
+Notice that $u(x) = u^h(x)$ at $x=0$ and $x=1$. Also note that
+\begin{equation}
+u_{,x}( \frac{1}{2} ) = u^h_{,x}( \frac{1}{2} ).
+\end{equation}
+
+\noindent {\bf Case 3:} $f(x) = p \, x ~$ where $p$ is a constant. In this
+case
+\begin{eqnarray}
+u(x) &=& g ~+~ (1 - x) \, h ~+~ \frac{p(1-x^3)}{6}\\
+u^h (x) \, &=& (h + g) \, ( 1 ~-~ x ) ~+~ g \, x + \frac{p(1-x)}{6}\nonumber \\
+ &=& g ~+~ (1 - x) \, h + \frac{p(1-x)}{6}
+\end{eqnarray}
+Notice that $u(x) = u^h(x)$ at $x=0$ and $x=1$. Also note that
+\begin{equation}
+u_{,x}( \frac{1}{\sqrt{3}} ) = u^h_{,x}( \frac{1}{\sqrt{3}} ).
+\end{equation}
+There are four points worth remembering from this exercise
+\begin{enumerate}
+\item The homogeneous equation ($f = 0$) was exact.
+
+\item With $f \ne 0$, $u^h = u$ at $x=0$ and $x=1$.
+
+\item With $f \ne 0$, there was at least one point where $u_{,x^h} = u{,_x}$
+
+\item The highest order polynomial in $x$ in the finite element solution $u^h$ is
+never greater than the highest order polynomial in the shape function.
+
+\end{enumerate}
+
+\subsection{Shape Functions}
+
+At this point, I will narrow the focus to deal with specifically the elements
+in my code, ConMan. It is possible to think very general shape functions,
+but in practice, people use triangles or quadralaterals (in 2-D). We will
+extend our finite element formulation to a 2-D equation next time.
+In terms of the level of approximation, there are also a lot of possibilities.
+We will stick to the simplest form, bi-linear elements, but you should be aware
+that higher order elements (bi-quadratic or bi-cubic spline elements) are also
+popular with some people. I am condensing a lot of very useful material from
+Chapter 3 of Hughes' book into one lecture. If you want to see more complete
+derivations, proofs of convergence, etc., of how to go about using higher order
+elements, look at Hughes book, Chapter 3.
+
+Lets start by thinking of a rectangle that is 2a by 2b in length centered at
+(0,0). There are two properties we would like the shape functions to have
+\begin{eqnarray}
+\sum_{A=1} ^4 N_A(X,Y) & = & 1 \label{eq:normalize} \\
+\sum_{A=1} ^4 N_A(X,Y) \, X_A & = & X \label{eq:interpx} \\
+\sum_{A=1} ^4 N_A(X,Y) \, Y_A & = & Y \label{eq:interpy}
+\end{eqnarray}
+Equation~\ref{eq:normalize} says that they are normalized, so that they sum to
+one (everywhere on X,Y). Equations ~\ref{eq:interpx} and~\ref{eq:interpy}
+state that the shape functions are also interpolation functions. Without
+doing a lot of derivation, I will claim that for the rectangle described
+above,
+\begin{eqnarray}
+N_1 & = & \frac{(a-x)(b-y)}{4ab} \\
+N_2 & = & \frac{(a+x)(b-y)}{4ab} \\
+N_3 & = & \frac{(a+x)(b+y)}{4ab} \\
+N_4 & = & \frac{(a-x)(b+y)}{4ab}
+\end{eqnarray}
+these shape functions satisfy the conditions in Equation~\ref{eq:normalize} and
+Equations ~\ref{eq:interpx} and~\ref{eq:interpy}. A good exercise would be to
+show this is true.
+
+\begin{quote}
+{\bf Problem 7} Show that the shape functions defined
+above satisfy Equation~\ref{eq:normalize} and Equations
+~\ref{eq:interpx} and~\ref{eq:interpy}.
+\end{quote}
+
+Notice that by convention, I start numbering my element nodes in the lower left
+hand corner and work counter-clockwise. {\em This is a very important point.
+It is followed through out all my finite element codes.} There is no magic
+reason, you just have to choose a starting place.
+
+In ConMan, as in Hughes, we further choose to normalize this by setting $a = 1$
+and $b=1$. This choice gives us an element whose area is 1, which is a
+convenient way to think about things. (This is because we left the factor of
+1/4 in the denominator). In my code, to make one less set of computations, I in
+effect set $a=0.5$ and $b=0.5$ so that the denominator goes to 1.
+
+Notice it is pretty easy to take derivatives of these shape functions
+\begin{eqnarray}
+N_{1,x} & = & \frac{-(1-y)}{4} \\
+N_{2,x} & = & \frac{(1-y)}{4} \\
+N_{3,x} & = & \frac{(1+y)}{4} \\
+N_{4,x} & = & \frac{-(1+y)}{4}
+\end{eqnarray}
+\begin{eqnarray}
+N_{1,y} & = & \frac{-(1-x)}{4} \\
+N_{2,y} & = & \frac{-(1+x)}{4} \\
+N_{3,y} & = & \frac{(1+x)}{4} \\
+N_{4,y} & = & \frac{(1-x)}{4}
+\end{eqnarray}
+What do we do if we want to solve a problem on a domain that is not convenient
+to split into a grid of 1 by 1 unit elements? We use an important principle
+of mathematics, the jacobian of the transformation
+\begin{equation}
+K_{11} ~=~ \int _A ^B N_{1,x} N_{1,x} \, dx ~=~ \int _0 ^1 N_{1,x} N_{1,x} \,
+J dx
+\end{equation}
+where $J$ is the Jacobian of the transformation. This is a very powerful
+point. When we are thinking of solving a regular Cartesian domain, this just
+corresponds to a stretching or a shrinking (notice we set $a=b=1$ above.
+However, if we are thinking about a cylindrical geometry, for example, we can
+use the Jacobian of the transformation between the geometries. Lets look at two
+examples:
+
+Converting an element 0.05 by 0.10 centered at (0.1,0.2) to the `parent element'
+centered at (0,0). Hughes also uses $\xi, \eta$ for the $X,Y$ coordinate
+pair in the `parent element' So we could write
+\begin{eqnarray}
+x & = & 0.1 + {0.05} \xi + 0.0 \eta\\
+y & = & 0.2 + {0.0} \xi + {0.10} \eta
+\end{eqnarray}
+or in matrix form we could write
+\begin{equation}
+\left \{ \begin{array}{c} x \\ y \end{array} \right \} =
+\left [ \begin{array}{cc}
+0.05 & 0.0 \\
+0.0 & 0.10
+\end{array} \right ] \left \{ \begin{array}{c}\xi \\ \eta \end{array} \right \}
++ \left \{ \begin{array}{c} 0.1 \\ 0.2 \end{array} \right \}
+\end{equation}
+Where $[J]$ is the Jacobian of the transformation. If the transformation were
+from an arbitratry shaped quadralateral to the parent element, then the off
+diagonal terms will not be zero. It is easy enough to show that
+\begin{equation}
+\int _{x_1} ^{x_2} \int _{y_1} ^{y_2} f(x,y) \, dx\, dy =
+\int _{-1} ^1 \int _{-1} ^1 f( \xi , \eta ) \, {\rm det} [J] \, d\xi \, d\eta
+\end{equation}
+It turns out, and it is also easy to show, that ${\rm det} [J]$ is the ratio of
+the areas when going from one rectangle to another (in fact any Cartesian to
+Cartesian transformation).
+
+{\bf Advanced Topic:} Now suppose we want to map a cylindrical domain to our
+`parent element.' We can use the same principle in this case:
+\begin{eqnarray}
+x &=& r \cos \theta = \cos \theta \, \xi - r \sin \theta \, \eta \\
+y &=& r \sin \theta = \sin \theta \, \xi + r \cos \theta \, \eta
+\end{eqnarray}
+so
+\begin{equation}
+{\rm det} [J_{geometry}] ~=~ r \cos^2 \theta + r \sin^2 \theta = r.
+\end{equation}
+of we would get
+\begin{equation}
+\int _{r_1} ^{r_2} \int _{\theta_1} ^{\theta_2} f(r\cos\theta, r\sin\theta) \,
+r \, dr\, d\theta =
+\int _{-1} ^1 \int _{-1} ^1 f( \xi , \eta ) \, {\rm det} [J_{area}] \, d\xi \,
+d\eta
+\end{equation}
+
+\subsubsection{Gauss Quadrature}
+
+An amazing fact, that makes the idea of finite elements easy and powerful is
+Gauss Quadrature. Gauss Quadrature is a way to turn an integral into a
+summation. Lets begin with a 1-D example, f(x) = c
+\begin{equation}
+\int _{-1} ^1 c \, dx = c x | _{-1} ^1 = 2 c
+\end{equation}
+Gauss noted that for any linear function
+\begin{equation}
+\int _{-1} ^{1} f(x) \, dx = 2.0 \times f( 0 ) = 2 c
+\end{equation}
+For a linear function, f(x) = a x + b
+\begin{equation}
+\int _{-1} ^{1} f(x) \, dx = f(\frac{-1}{\sqrt{3}}) + f(\frac{-1}{\sqrt{3}})
+\end{equation}
+The direct way,
+\begin{equation}
+\int _{-1} ^{1} (ax + b) \, dx = (\frac{ax^2}{2} + bx) | _{-1} ^1 =
+\frac{a}{2} + b - (\frac{a}{2} -b) = 2b
+\end{equation}
+Gauss' way
+\begin{equation}
+\int _{-1} ^{1} (ax + b) \, dx = a \frac{-1}{\sqrt{3}} + b + a \frac{1}{\sqrt{3}}
++ b = 2b
+\end{equation}
+It turns out that $\frac{-1}{\sqrt{3}}, \frac{1}{\sqrt{3}}$ are exact for a
+linear equation, but from what I showed, so would any $-x, x$ combination, but
+what Gauss showed was more powerful, that if the function is of higher order,
+the $ \frac{-1}{\sqrt{3}}, \frac{1}{\sqrt{3}}$ choice is the best approximation
+you can make with only two terms. If we go to three terms, the choice would
+be $ -\sqrt{\frac{3}{5}} , 0 , \sqrt{\frac{3}{5}}$.
+
+To integrate a 2-D Cartesian region, like our parent element, it turns out that
+2 by 2 quadrature, or the four points
+\begin{eqnarray}
+\xi = \frac{-1}{\sqrt{3}} & ~~~ & \eta = \frac{-1}{\sqrt{3}} \\
+\xi = \frac{ 1}{\sqrt{3}} & ~~~ & \eta = \frac{-1}{\sqrt{3}} \\
+\xi = \frac{ 1}{\sqrt{3}} & ~~~ & \eta = \frac{ 1}{\sqrt{3}} \\
+\xi = \frac{-1}{\sqrt{3}} & ~~~ & \eta = \frac{ 1}{\sqrt{3}}
+\end{eqnarray}
+are sufficient to exactly integrate our bilinear shape functions over the
+{-1,-1} to {1,1} domain.
+
+At this point, it would be worth talking about the code ConMan for a minute.
+The shape functions are generated in ConMan in the routine {\bf genshp} for
+GENerate SHape functions Parent domain. If you look at the routine you will
+find the first part of it is pretty easy to follow from the discussion in
+todays lecture. Some of the second part is a little tricky in its details,
+but generally it is also pretty easy to follow. The subroutine local will
+be discussed in detail next lecture. The subroutine genshg will follow this
+one. Some constants to keep in mind (ConMan was written very generally,
+so we tried not to hard-wire element sizes into it.
+\begin{enumerate}
+\item numnp = total NUMber of Nodal Points
+
+\item numel = total NUMber of ELements
+
+\item nsd = Number of Spacial Dimensions (2 for a 2-D problem).
+
+\item nen = Number of Element Nodes (4 for the bi-linear quad).
+
+\item shl = SHape function Local (array corresponding to equations 54-57).
+
+\item shdx = SHape function Derivative X (array of x derivs)
+
+\item shdy = SHape function Derivative Y (array of y derivs)
+\end{enumerate}
+
+\tt
+\begin{verbatim}
+ subroutine genshp ( x, xl ,ien , locblk , shl , shdx ,
+ & shdy, det, nel , nblk , elval)
+c
+c----------------------------------------------------------------------
+c
+c This program generates the shape functions for bi-linear,and
+c calculates the min element dimension per node
+c
+c input:
+c x (nsd,numnp) : coordinates
+c xl (lvec,nsd,nen) : local coordinates
+c ien (numel,nen) : ien array
+c nel : number of elements
+c nblk : number of element blocks
+c
+c output:
+c
+c shl (nen,nipt)
+c shdx (nel,nen,nipt)
+c shdy (nel,nen,nipt)
+c det (nel,nipt)
+c eval (nel,6)
+c
+c Note: the last five arrays are setup with element as the first index.
+c
+c
+ implicit double precision (a-h,o-z)
+c
+c.... deactivate above card(s) for single precision operation
+c
+ include 'common.h'
+c
+ dimension x(nsd,1) , ien(nel), locblk(2,1),
+ & shl(nen,1) , shdx(1) , shdy(1) ,
+ & det(1) , xl(lvec,nsd,nen) , elval(nel,1)
+c
+ common /temp1 / sa(4),ta(4),sg(5),tg(5),shldx(4,5),shldy(4,5)
+c
+c.... set up parameters
+c
+ sa(1) = -pt5
+ sa(2) = pt5
+ sa(3) = pt5
+ sa(4) = -pt5
+c
+ ta(1) = -pt5
+ ta(2) = -pt5
+ ta(3) = pt5
+ ta(4) = pt5
+c
+ guass = one / sqrt(three)
+c
+
+ sg(1) = -one * guass
+ sg(2) = one * guass
+ sg(3) = one * guass
+ sg(4) = -one * guass
+ sg(5) = zero * guass
+c
+ tg(1) = -one * guass
+ tg(2) = -one * guass
+ tg(3) = one * guass
+ tg(4) = one * guass
+ tg(5) = zero * guass
+c
+ call clear (shdx, nel*nen*nipt)
+ call clear (shdy, nel*nen*nipt)
+c
+c generate the generic shape functions
+c
+ do 100 j = 1,5
+ do 100 i = 1,4
+ shl(i,j) = (pt5 + sa(i) * sg(j) ) * (pt5 + ta(i) * tg(j) )
+ shldx(i,j) = sa(i)*(pt5 + ta(i) * tg(j) )
+ shldy(i,j) = ta(i)*(pt5 + sa(i) * sg(j) )
+100 continue
+c
+ do 1000 iblk = 1, nblk
+c
+c.... set up the pointers
+c
+ iel = locblk(1,iblk)
+ nenl = locblk(2,iblk)
+ nvec = locblk(1,iblk+1) - iel
+ call local (ien(iel), x, xl, nsd, nenl, 'localize')
+c
+c.... call the global shape function
+c
+ call genshg (shl , shldx , shldy , shdx,
+ & shdy , det , xl , iel )
+c
+c.... calculate min element dimension per node
+c
+ do 500 iv = 1 ,nvec
+ ivel = iel + iv - 1
+ exse1 = pt5*( xl(iv,1,2) + xl(iv,1,3) - xl(iv,1,4) - xl(iv,1,1) )
+ exse2 = pt5*( xl(iv,2,2) + xl(iv,2,3) - xl(iv,2,4) - xl(iv,2,1) )
+ eeta1 = pt5*( xl(iv,1,3) + xl(iv,1,4) - xl(iv,1,1) - xl(iv,1,2) )
+ eeta2 = pt5*( xl(iv,2,3) + xl(iv,2,4) - xl(iv,2,1) - xl(iv,2,2) )
+ hxse = sqrt(exse1*exse1 + exse2*exse2)
+ heta = sqrt(eeta1*eeta1 + eeta2*eeta2)
+ elval(ivel ,1) = exse1/hxse
+ elval(ivel ,2) = exse2/hxse
+ elval(ivel ,3) = eeta1/heta
+ elval(ivel ,4) = eeta2/heta
+ elval(ivel ,5) = hxse
+ elval(ivel ,6) = heta
+500 continue
+c
+1000 continue
+c
+c.... return
+c
+ return
+ end
+\end{verbatim}
+
+\rm
+Next we look at {\bf genshg} for GENerate SHape functions Global. This is
+where the parent elements get mapped to the `real' elements. Basically it
+calculates a lot of determinents.
+
+\tt
+\begin{verbatim}
+ subroutine genshg (shl , shldx , shldy , shdx , shdy ,
+ & det , xl , iel )
+c
+c----------------------------------------------------------------------
+c
+c This program generates the globalshape functions for bi-linear,
+c
+c input:
+c xl (lvec,nsd,nen) : local coordinates
+c shldx (nen ,nipt ) : local dx
+c shldy (nen ,nipt ) : local dy
+c shl (nen ,nipt ) : local shape functions
+c c output:
+c
+c shdx (nel,nen,nipt)
+c shdy (nel,nen,nipt)
+c det (nel,nipt)
+c
+c Note: the last four arrays are setup with element as the first index.
+c This should facilitate vectorization.
+c
+c
+c----------------------------------------------------------------------
+c
+c
+ implicit double precision (a-h,o-z)
+c
+c.... deactivate above card(s) for single precision operation
+c
+ include 'common.h'
+c
+ dimension shl(nen,1), shldx(nen,1), shldy(nen,1), xl(lvec,nsd,1)
+ dimension shdx(numel,nen,1) ,shdy(numel,nen,1) ,det(numel,1)
+c
+ common /tempx/ xs(lvec,2,2) , temp(lvec)
+c
+c.... loop over all the integration points
+c
+ do 1000 l = 1 , nipt
+c
+c find jocabian
+c
+ do 100 iv = 1 , nvec
+ xs(iv,1,1) = xl(iv,1,1) * shldx(1,l) + xl(iv,1,2) * shldx(2,l)
+ & + xl(iv,1,3) * shldx(3,l) + xl(iv,1,4) * shldx(4,l)
+c
+ xs(iv,1,2) = xl(iv,2,1) * shldx(1,l) + xl(iv,2,2) * shldx(2,l)
+ & + xl(iv,2,3) * shldx(3,l) + xl(iv,2,4) * shldx(4,l)
+c
+ xs(iv,2,1) = xl(iv,1,1) * shldy(1,l) + xl(iv,1,2) * shldy(2,l)
+ & + xl(iv,1,3) * shldy(3,l) + xl(iv,1,4) * shldy(4,l)
+c
+ xs(iv,2,2) = xl(iv,2,1) * shldy(1,l) + xl(iv,2,2) * shldy(2,l)
+ & + xl(iv,2,3) * shldy(3,l) + xl(iv,2,4) * shldy(4,l)
+100 continue
+c
+c..... calculate the inverse jacobian
+c
+ do 200 iv = 1 , nvec
+ ivel = iv + iel - 1
+ det(ivel,l) = xs(iv,1,1) * xs(iv,2,2) - xs(iv,1,2) * xs(iv,2,1)
+200 continue
+c
+c.... check for zero determine
+c
+ do 300 iv = 1 , nvec
+ ivel=iv+iel-1
+ if ( det(ivel,l) .le. zero )
+ & call error (' genshg ','det-jacb',iel)
+300 continue
+c
+c.... continue inverse calculation
+c.... and find derivative with respect to global axes
+c
+ do 400 iv = 1 , nvec
+ ivel=iv+iel-1
+ temp(iv) = det(ivel,l)
+ temp(iv) = one/temp(iv)
+ shdx(ivel,1,l) = temp(iv) * ( xs(iv,2,2) * shldx(1,l)
+ & - xs(iv,1,2) * shldy(1,l) )
+ shdy(ivel,1,l) = temp(iv) * (-xs(iv,2,1) * shldx(1,l)
+ & + xs(iv,1,1) * shldy(1,l) )
+c
+ shdx(ivel,2,l) = temp(iv) * ( xs(iv,2,2) * shldx(2,l)
+ & - xs(iv,1,2) * shldy(2,l) )
+ shdy(ivel,2,l) = temp(iv) * (-xs(iv,2,1) * shldx(2,l)
+ & + xs(iv,1,1) * shldy(2,l) )
+c
+ shdx(ivel,3,l) = temp(iv) * ( xs(iv,2,2) * shldx(3,l)
+ & - xs(iv,1,2) * shldy(3,l) )
+ shdy(ivel,3,l) = temp(iv) * (-xs(iv,2,1) * shldx(3,l)
+ & + xs(iv,1,1) * shldy(3,l) )
+c
+ shdx(ivel,4,l) = temp(iv) * ( xs(iv,2,2) * shldx(4,l)
+ & - xs(iv,1,2) * shldy(4,l) )
+ shdy(ivel,4,l) = temp(iv) * (-xs(iv,2,1) * shldx(4,l)
+ & + xs(iv,1,1) * shldy(4,l) )
+400 continue
+c
+c.... end of integration point loop
+c
+1000 continue
+c
+c.... return
+c
+ return
+ end
+\end{verbatim}
+\rm
+\clearpage
+
+There are two strange data structures and loops that appear in the subroutines
+above. The first is,
+\tt\begin{verbatim}
+ do 1000 iblk = 1, nblk
+c
+c.... set up the pointers
+c
+ iel = locblk(1,iblk)
+ nenl = locblk(2,iblk)
+ nvec = locblk(1,iblk+1) - iel
+\rm
+\end{verbatim}
+and the other is loops like
+\tt\begin{verbatim}
+ do 200 iv = 1 , nvec
+ ivel = iv + iel - 1
+ det(ivel,l) = xs(iv,1,1) * xs(iv,2,2) - xs(iv,1,2) * xs(iv,2,1)
+200 continue
+\end{verbatim}\rm
+These are done because on a vector computer you can not be sure that assemble
+operations will be handled correctly unless you are very careful. This
+involves a little thinking into how a vector computer works. I have removed all of this
+coding from the test code you will be working with, so don't get hung up on it.
+It is a relic from the days of Cray computers.
+
+\clearpage
+
+\subsection{The Element Point of View}
+
+Transforming between the element and global points of view is done
+with the data structure called the IEN array for Element to Node
+transformation. The IEN array takes an element number and a local node
+number and it's value is the global node number. It is easiest to look at
+some examples:
+
+Consider a 3 element by 3 element grid. I will number my elements and node
+starting in the lower left hand corner, and working up fastest.
+
+\begin{verbatim}
+For element 3:
+ ien (3, 1) = 3
+ ien (3, 2) = 7
+ ien (3, 3) = 8
+ ien (3, 4) = 4
+For element 5:
+ ien (5, 1) = 6
+ ien (5, 2) = 10
+ ien (5, 3) = 11
+ ien (5, 4) = 7
+\end{verbatim}
+
+There are three kinds of opperations we might think of wanting. The first is
+taking values of some function (coordinates, velocities, temperatures,
+stresses, etc.) defined on the global grid and getting the values for a single
+element, this is called a {\em gather} operation. The next is taking values in
+an element and spreading them out to the global array, this is called a {\em
+scatter} operation. The third operation is to take the value at a local node
+and add it to the global value for that node, an {\em assembly} step.
+
+All three operations, gather, scatter, and assemble are done by the routine
+local. Because it is called by the genshp routine above, it is a good example
+to look at. Next class time I will show you how it looks in the code where we
+generate the stiffness matrix.
+
+\tt
+\begin{verbatim}
+ subroutine local (ien, global, rlocal, n, nenloc, code)
+c
+c----------------------------------------------------------------------
+c
+c This subroutine localizes the values via the element group IEN
+c array. The localize values are spread with the element number as
+c its first (fastest) index for vectorization.
+c
+c input:
+c ien (numel,nen) : the element connectivity
+c global (n,numnp) : the global array
+c rlocal (lvec,n,nenl) : the local array
+c n : number of d.o.f's to be copied
+c nenl : local number of element nodes
+c code : the transfer code
+c .eq. 'localize', from global to local
+c .eq. 'globaliz', from local to global
+c .eq. 'add glob', add local to global
+c nvec : vector length, number of elements ( pasted from common)
+c nelvec : number of elments in group ( pasted from common)
+c
+c output:
+c----------------------------------------------------------------------
+c
+c
+ implicit double precision (a-h,o-z)
+c
+c.... deactivate above card(s) for single precision operation
+c
+ include 'common.h'
+c
+ dimension rlocal(lvec,n,1), global(n,1), ien(numel,1)
+ character*8 code
+c
+c.... --------------------> 'localize' <--------------------
+c
+ if (code .eq. 'localize') then
+c
+c.... loop over element nodes, d.o.f's and number of element
+c
+c$dir scalar
+ do 300 ienl = 1, nenloc
+c$dir scalar
+ do 200 i = 1, n
+c$dir no_recurrence
+ do 100 iv = 1, nvec
+ rlocal(iv,i,ienl) = global(i,ien(iv,ienl))
+100 continue
+200 continue
+300 continue
+c
+c.... return
+c
+ return
+ endif
+c
+c.... --------------------> 'globaliz' <--------------------
+c
+ if (code .eq. 'globaliz') then
+c
+c.... loop over element nodes, d.o.f's and number of element
+c
+c$dir scalar
+ do 600 ienl = 1, nenloc
+c$dir scalar
+ do 500 i = 1, n
+c$dir no_recurrence
+ do 400 iv = 1, nvec
+ global(i,ien(iv,ienl)) = rlocal(iv,i,ienl)
+400 continue
+500 continue
+600 continue
+c
+c.... return
+c
+ return
+ endif
+c
+c.... --------------------> 'add glob' <--------------------
+c
+ if (code .eq. 'add glob') then
+c
+c.... loop over element nodes, d.o.f's and number of element
+c
+c$dir scalar
+ do 900 ienl = 1, nenloc
+c$dir scalar
+ do 800 i = 1, n
+c$dir no_recurrence
+ do 700 iv = 1, nvec
+ global(i,ien(iv,ienl)) = global(i,ien(iv,ienl)) +
+ & rlocal(iv,i,ienl)
+700 continue
+800 continue
+900 continue
+c
+c.... return
+c
+ return
+ endif
+c
+c.... --------------------> error <--------------------
+c
+ call error ('local ', code, 0)
+c
+c.... end
+c
+ end
+\end{verbatim}
+\rm
+
+In practice, one would like to separate the gather/scatter and assemble
+opperations from the bulk of the computation. One very important reason for
+this is because some computers (vector computers like Crays and parallel
+computers) require special thought in how one gathers and scatters information
+(because more than one processor may be working on part of the problem at the
+same time). I don't want to add further confusion at this point by diving into
+that problem, but that is what ConMan is doing at certain points. In addition,
+you will find as you go through ConMan that I seldom use the routine local.
+Often I put the operations from local directly into the code. This is because
+at one time, it was popular to take any small segment of code and make it a
+subroutine. As the speed of floating point operations (multiplies and adds)
+became faster, the overhead associated with calling a subroutine bacame
+increasingly expensive compared to the cost of a floating point operation. Thus
+code writing moved away from having many small subrouines (as is the style of
+{\bf dlearn,} the code in Hughes' book, and more like ConMan.
+
+\subsection{Assembly}
+
+Boundary conditions and equation numbers (the lm array).
+
+\section{Stiffness Matrix for Viscous Flow: Penalty Method}
+
+There are several methods for solving viscous flow problems with the finite
+element method. 'Mixed methods' solve for velocities and pressures. They
+have the drawback that the amount of computation and the amount of storage
+increases. The 'penalty method' is a classical method that approaches the
+incompressible limit. A number of tests with the penalty formulation have
+shown that the solutions are in very good agreement with solutions from other
+methods, such as stream-function methods which are truely incompressible. The
+penalty method eliminates the pressure from the equations, but there are some
+minor problems that result (nothing comes for free). If you want to look at
+the results of some 2-D incompressible codes and see how the penalty method
+compares, you should look at the following references.
+
+\begin{quotation}
+\noindent
+Travis, B. J., C. Anderson, J. Baumgardner, C. Gable, B. H. Hager, P.
+Olson, R. J. O'Connell, A. Raefsky, and G. Schubert, A benchmark comparison of
+numerical methods for infinite Prandtl number convection in two-dimensional
+Cartesian geometry. {\em Geophys. Astrophys. Fluid Dyn., 55,} 137-180, 1990.
+
+\medskip
+\noindent
+Blankenbach, B., et al., A benchmark comparison for mantle convection codes,
+{\it Geophys. J. Int., 98,} 23-38, 1989.
+\end{quotation}
+
+At this point, lets return to the 2-D viscous flow equations, written in terms
+of stresses, instead of velocities
+\begin{eqnarray}
+\tau_{ij,j} + f_i & = & 0 \label{eq:motion}\\
+u_{i,i} & = & 0
+\end{eqnarray}
+where
+\begin{equation}
+\tau_{ij} = -p \delta_{ij} + 2 \mu u_{(i,j)} \label{eq:constit}
+\end{equation}
+where
+\begin{equation}
+u_{(i,j)} = (u_{i,j} + u_{j,i})/2
+\end{equation}
+We replace equation~\ref{eq:constit} with the following relationships
+\begin{eqnarray}
+\tau_{ij} = -p^\lambda \delta_{ij} + 2 \mu u_{(i,j)} \label{eq:constit-lam} \\
+0 = u_{i,i} + p^\lambda/\lambda. \label{eq:constit-lam2}
+\end{eqnarray}
+As $\lambda$ approaches infinity, these relations approach the incompressible
+solution. Also, as $\lambda$ approaches infinity, $p^\lambda$ approaches the
+hydrostatic pressure in the incompressible case. In general, the hydrostatic
+pressure is $-\tau_{ii}/3$. Substituting Equation~\ref{eq:constit-lam2}
+into~\ref{eq:constit-lam} we get
+\begin{equation}
+\tau_{ij} = \lambda u_{i,i} \delta_{ij} + 2 \mu u_{(i,j)}
+\label{eq:constit-lam-final}
+\end{equation}
+or
+\begin{equation}
+\tau_{ii} = 3 \lambda u_{i,i} + 2 \mu u_{i,i}
+\end{equation}
+or
+\begin{equation}
+\tau_{ii}/3 = -p = (\lambda + 2/3 \mu) u_{i,i}
+\end{equation}
+but we also have
+\begin{equation}
+-p^\lambda = \lambda u_{i,i}
+\end{equation}
+from Equation~\ref{eq:constit-lam2}. Clearly in the incompressible limit
+$\lambda \gg \mu$ then $\lambda + 2/3 \mu \rightarrow \lambda$ and $p^\lambda
+\rightarrow p$. Also note that the continuity equation is satisfied.
+
+Now, substituting Equation~\ref{eq:constit-lam-final} into
+Equation~\ref{eq:motion} we have
+\begin{equation}
+\{ \lambda u_{i,i} \delta_{ij} + 2 \mu u_{(i,j)} \},j + f_i = 0
+\end{equation}
+At this point, it is probably easier to switch to differential notation. I
+will also specialize to 2-D:
+\begin{eqnarray}
+\frac{\partial}{\partial x} \{
+\lambda ( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial z})
++ 2 \mu (\frac{\partial u}{\partial x} + \frac{\partial u}{\partial x} )/2 \}+
+\frac{\partial}{\partial z} \{
+2 \mu (\frac{\partial v}{\partial x} + \frac{\partial u}{\partial z} )/2 \}+
+f_x = 0 \\
+\frac{\partial}{\partial x} \{
+2 \mu (\frac{\partial u}{\partial z} + \frac{\partial v}{\partial x} )/2 \}+
+\frac{\partial}{\partial z} \{
+\lambda ( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial z} )
++ 2 \mu (\frac{\partial v}{\partial z} + \frac{\partial v}{\partial z} )/2 \}+
+f_z = 0
+\end{eqnarray}
+
+These are second order partial differential equations. Simplifying, I get
+\begin{eqnarray}
+\lambda (
+\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 v}{\partial x\partial z} )
++ 2 \mu \frac{\partial^2 u}{\partial x^2}
++ \mu ( \frac{\partial^2 u}{\partial z^2}
++ \frac{\partial^2 v}{\partial z\partial x} )
++ f_x = 0 \\
+\lambda (
+\frac{\partial^2 u}{\partial z\partial x} + \frac{\partial^2 v}{\partial z^2} )
++ \mu ( \frac{\partial^2 u}{\partial x\partial z}
++ \frac{\partial^2 v}{\partial x^2} )
++ 2 \mu \frac{\partial^2 v}{\partial z^2}
++ f_z = 0
+\end{eqnarray}
+
+Now we use the same technique (approach) as we used in Possion's equation to
+turn the differential form into an integral form. You can either look at it
+as we find the variational form of the stokes equation (which is what we are
+doing) or you can think of it as multiplying by a weighting function $w$ and
+integrating over the domain. Then using integration by parts to convert the
+second derivatives to first derivatives. This is done in carefully by
+Hughes on pages 197-200, but he has left out a number of intermediate steps.
+Nothing about this step is hard, it is just tedious. There is, however, a
+cleaver short cut. If we return to the messy equations at the top of the
+page, multiply them by the weighting function $w$ and integrate over the
+domain, then we do not have to use integration by parts. If you are confused,
+or don't believe me, then you should take the equations directly above this
+paragraph, multiply by a weighting function $w$ and integrate over the 2-D
+domain $\Omega$, then use integration by parts. You will find (after a little
+algebra)
+
+\begin{eqnarray}
+\int\int _{\Omega} \frac{\partial w}{\partial x} \{
+\lambda ( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial z})
++ 2 \mu \frac{\partial u}{\partial x} \}+
+\frac{\partial w}{\partial z} \{ 2 \mu (\frac{\partial v}{\partial x} +
+\frac{\partial u}{\partial z} )/2 \} \, d\Omega + \nonumber \\
+\int\int _{\Omega} f_x \, w \, d\Omega = b.c.~terms \\
+\int\int _{\Omega} \frac{\partial w}{\partial x} \{
+2 \mu (\frac{\partial u}{\partial z} + \frac{\partial v}{\partial x} )/2 \}+
+\frac{\partial w}{\partial z} \{
+\lambda ( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial z} ) +
+2 \mu \frac{\partial v}{\partial z} \} \, d\Omega + \nonumber \\
+\int\int _{\Omega} f_z \, w \, d\Omega = b.c.~terms
+\end{eqnarray}
+Note that we don't get something for nothing, this short cut does not give us
+the boundary condition terms (velocity or flux). These would fall out of the
+integration by parts. Recall,
+\begin{equation}
+\int _a ^b w \, dv = w\, v | _a ^b - \int _a ^b v \, dw
+\end{equation}
+were in our case $w$ is the weighting function and $v$ is the second derivative
+term. The first term gives us the flux (first derivative) boundary
+conditions. In the case of the momentum equations, that is the applied
+tractions (or stress boundary conditions).
+
+Now we make use of Galerkin's approximation, or more simply, we use the same
+weighting functions as we use for interpolation function, i.e., the shape
+functions, N. So we substitute
+\begin{eqnarray}
+\frac{\partial w}{\partial x} = N_x \\
+\frac{\partial w}{\partial z} = N_z \\
+\frac{\partial u}{\partial x} = u\, N_x \\
+\frac{\partial u}{\partial z} = u\, N_z \\
+\frac{\partial v}{\partial x} = v\, N_x \\
+\frac{\partial v}{\partial z} = v\, N_z
+\end{eqnarray}
+into our weak form equations. Although messy, that is straight-forward.
+\begin{eqnarray}
+\int\int _{\Omega} N_x \{
+\lambda ( u\, N_x + v\, N_z) + 2 \mu u\, N_x \} +
+N_z \{ \mu ( v\, N_x + u\, N_z ) \} \, d\Omega + \nonumber \\
+\int\int _{\Omega} f_x \, w \, d\Omega = b.c.~terms \\
+\int\int _{\Omega} N_x \{ \mu ( u\, N_z + v\, N_x ) \}+
+N_z \{ \lambda ( u\, N_x + v\, N_z ) +
+2 \mu v\, N_z \} \, d\Omega + \nonumber \\
+\int\int _{\Omega} f_z \, w \, d\Omega = b.c.~terms
+\end{eqnarray}
+
+At this point, it is useful to separate the equations into a $\lambda$ part and
+a $\mu$ part. We can also write them as a 2-D matrix equation
+\begin{equation} [K_{\lambda}] = \left[
+\begin{array}{cc}
+N_x \lambda N_x & N_x \lambda N_z \\
+N_z \lambda N_x & N_z \lambda N_z
+\end{array} \right ] \label{eq:Klambda}
+\end{equation}
+and
+\begin{equation}
+[K_{\mu}] = \left[
+\begin{array}{cc}
+N_x 2 \mu N_x + N_z \mu N_z & N_z \mu N_x \\
+N_x \mu N_z & N_z 2 \mu N_z + N_x \mu N_x
+\end{array} \right ]. \label{eq:Kmu}
+\end{equation}
+
+Hughes makes use of an interesting, and important observation. This
+observation will greatly simplify constructing the stiffness matrix for
+arbitrary coordinate systems. We can rewrite the stiffness matrices above in
+the following form:
+\begin{equation}
+[K_{\lambda}] + [K_{\mu}] = [B]^T [D] [B]
+\end{equation}
+\begin{equation}
+[D_{\lambda}] + [D_{\mu}] = [D]
+\end{equation}
+where
+\begin{equation}
+[D_{\mu}] = \mu \left[
+\begin{array}{ccc}
+2 & 0 & 0 \\
+0 & 2 & 0 \\
+0 & 0 & 1
+\end{array} \right ]
+\end{equation}
+and
+\begin{equation}
+[D_{\lambda}] = \lambda \left[
+\begin{array}{ccc}
+1 & 1 & 0 \\
+1 & 1 & 0 \\
+0 & 0 & 0
+\end{array} \right ]
+\end{equation}
+and
+\begin{equation}
+[B] = \left[
+\begin{array}{cc}
+N_x & 0 \\
+0 & N_z \\
+N_z & N_x
+\end{array} \right ].
+\end{equation}
+
+Notice what the $[B]$ matrix is...
+\begin{equation}
+[B] \{\vec u\} = ?
+\end{equation}
+
+\subsection{The Element Stiffness Matrix For 2-D Viscous Flow}
+
+Lets return to our 2-D, four node, bi-linear-velocity element. We want to
+form the `element stiffness matrix' for that element. Then we will assemble
+the element stiffness matrices for all the elements together to form the global
+stiffness matrix. Because there are 2 degrees of freedom ($u,v$) for each
+node, there will be eight unknowns for each element. In the same way that
+the 2 by 2 element matricies are symmetric, it will turn out that element
+stiffness matrix is symmetric, and infact the global matrix will be symmetric.
+Symmetric matricies are one of the properties of Galerkin formulations, you
+can read the details in Hughes. Because we want to be stingy with memory, we
+will only form and store the upper triangular part of the stiffness matrix.
+If we formed these one element at a time, the savings would not be worth the
+effort, but because we will form many elements at one time on a vector machine,
+we need to be careful. In ConMan, I treat the 8 by 8 stiffness matrix as a 1-D
+array with 36 elements. The numbering as as follows:
+\begin{equation}
+[K^e] = \left[
+\begin{array}{cccccccc}
+ 1 & 2 & 4 & 7 & 11 & 16 & 22 & 29 \\
+ & 3 & 5 & 8 & 12 & 17 & 23 & 30 \\
+ & & 6 & 9 & 13 & 18 & 24 & 31 \\
+ & & & 10 & 14 & 19 & 25 & 32 \\
+ & & & & 15 & 20 & 26 & 33 \\
+ & & & & & 21 & 27 & 34 \\
+ & & & & & & 28 & 35 \\
+ & & & & & & & 36
+\end{array} \right ].
+\end{equation}
+You can think of this 8 by 8 matrix as sixteen 2 by 2 matricies. Each 2 by 2
+matrix is made of elements from the 2 by 2 matrices~\ref{eq:Klambda}
+and~\ref{eq:Kmu}.
+
+Recall that the shape functions for the four nodes look like
+\begin{eqnarray}
+N_1 & = & \frac{(a-x)(b-y)}{4ab} \\
+N_2 & = & \frac{(a+x)(b-y)}{4ab} \\
+N_3 & = & \frac{(a+x)(b+y)}{4ab} \\
+N_4 & = & \frac{(a-x)(b+y)}{4ab}
+\end{eqnarray}
+or in the fortran code, these are stored in the array shl as: shl(1,gp),
+shl(2,gp), shl(3,gp), shl(4,gp), where {\em gp} is the number of the gauss
+point.
+
+If we go back to the 2 by 2 matrix for the $\lambda$ part of the stiffness
+matrix, Equation~\ref{eq:Kmu}, note that the first shape function in each
+matrix element came from the weighting function and the second shape function
+came from the interolation function
+\begin{equation} [K_{\lambda}] = \left[
+\begin{array}{cc}
+N_x(from~w) \lambda N_x(from~u_g) & N_x(from~w) \lambda N_z(from~v_g) \\
+N_z(from~w) \lambda N_x(from~u_g) & N_z(from~w) \lambda N_z(from~v_g)
+\end{array} \right ]
+\end{equation}
+This helps illustrate how the elements of the 8 by 8 matrix come together. I'm
+going to substitute things in one step at a time, to help you visualize what is
+happening. Lets substitue the fortran arrays in for the shape function
+derivatives for the interpolation functions, and the $u_g$'s.
+\begin{eqnarray}
+N_x(from~u_g) & = & shdx(ivel,i,gp) \\
+N_z(from~u_g) & = & shdy(ivel,i,gp) \\
+u_g & = & vl(iv,k), where~k=1,3,5,~or~7 \\
+v_g & = & vl(iv,k), where~k=2,4,6,~or~8
+\end{eqnarray}
+
+\begin{displaymath}
+[K^e_{\lambda}] =
+\end{displaymath}
+{\tiny
+\begin{displaymath}\left[
+\begin{array}{cccc}
+N_x \lambda shdx(ivel,1,gp) & N_x \lambda shdy(ivel,1,gp) &
+N_x \lambda shdx(ivel,2,gp) & N_x \lambda shdy(ivel,2,gp) \\
+
+\null & N_z \lambda shdy(ivel,1,gp) &
+N_z \lambda shdx(ivel,2,gp) & N_z \lambda shdy(ivel,2,gp) \\
+
+\null & \null &
+N_x \lambda shdx(ivel,2,gp) & N_x \lambda shdy(ivel,2,gp) \\
+
+\null & \null &
+\null & N_z \lambda shdy(ivel,2,gp) \\
+
+\null & \null &
+\null & \null \\
+
+\null & \null &
+\null & \null \\
+
+\null & \null &
+\null & \null \\
+
+\null & \null &
+\null & \null \\
+\end{array} \right .
+\end{displaymath}
+\begin{equation} \left .
+\begin{array}{cccc}
+N_x \lambda shdx(ivel,3,gp) & N_x \lambda shdy(ivel,3,gp) &
+N_x \lambda shdx(ivel,4,gp) & N_x \lambda shdy(ivel,4,gp) \\
+
+N_z \lambda shdx(ivel,3,gp) & N_z \lambda shdy(ivel,3,gp) &
+N_z \lambda shdx(ivel,4,gp) & N_z \lambda shdy(ivel,4,gp) \\
+
+N_x \lambda shdx(ivel,3,gp) & N_x \lambda shdy(ivel,3,gp) &
+N_x \lambda shdx(ivel,4,gp) & N_x \lambda shdy(ivel,4,gp) \\
+
+N_z \lambda shdx(ivel,3,gp) & N_z \lambda shdy(ivel,3,gp) &
+N_z \lambda shdx(ivel,4,gp) & N_z \lambda shdy(ivel,4,gp) \\
+
+N_x \lambda shdx(ivel,3,gp) & N_x \lambda shdy(ivel,3,gp) &
+N_x \lambda shdx(ivel,4,gp) & N_x \lambda shdy(ivel,4,gp) \\
+
+\null & N_z \lambda shdy(ivel,3,gp) &
+N_z \lambda shdx(ivel,4,gp) & N_z \lambda shdy(ivel,4,gp) \\
+
+\null & \null &
+N_x \lambda shdx(ivel,4,gp) & N_x \lambda shdy(ivel,4,gp) \\
+
+\null & \null &
+\null & N_z \lambda shdy(ivel,4,gp) \\
+\end{array} \right ].
+\end{equation} }
+
+Think of this as the $\lambda$ part of the the element matrix equation
+\begin{equation}
+[ K^e ] \{ vl \} = \{ f^e \}
+\end{equation}
+where $\{ vl \}$ and $\{ f^e \}$ are the arrays of local unknowns and local
+buoyancy forces. In the same way that we substituted the shape functions for
+the interpolation function, we substitute the shape functions for the weighting
+functions (the remaining $N_x, N_z$'s in the matrix above. I am doing this
+step by step so you can see which node gets connected with which node....
+
+\begin{displaymath} [K^e_{\lambda}] =
+\end{displaymath} {\tiny
+\begin{displaymath}\left[
+\begin{array}{cccc}
+
+shdx(ivel,1,gp) \lambda shdx(ivel,1,gp) & shdx(ivel,1,gp) \lambda
+shdy(ivel,1,gp) &
+shdx(ivel,1,gp) \lambda shdx(ivel,2,gp) & shdx(ivel,1,gp) \lambda
+shdy(ivel,2,gp) \\
+
+\null & shdy(ivel,1,gp) \lambda shdy(ivel,1,gp) &
+shdy(ivel,1,gp) \lambda shdx(ivel,2,gp) & shdy(ivel,1,gp) \lambda
+shdy(ivel,2,gp) \\
+
+\null & \null &
+shdx(ivel,2,gp) \lambda shdx(ivel,2,gp) & shdx(ivel,2,gp) \lambda
+shdy(ivel,2,gp) \\
+
+\null & \null &
+\null & shdy(ivel,2,gp) \lambda shdy(ivel,2,gp) \\
+
+\null & \null &
+\null & \null \\
+
+\null & \null &
+\null & \null \\
+
+\null & \null &
+\null & \null \\
+
+\null & \null &
+\null & \null \\
+\end{array} \right .
+\end{displaymath}
+\begin{equation} \left .
+\begin{array}{cccc}
+shdx(ivel,1,gp) \lambda shdx(ivel,3,gp) & shdx(ivel,1,gp) \lambda
+shdy(ivel,3,gp) &
+shdx(ivel,1,gp) \lambda shdx(ivel,4,gp) & shdx(ivel,1,gp) \lambda
+shdy(ivel,4,gp)
+\\
+
+shdy(ivel,1,gp) \lambda shdx(ivel,3,gp) & shdy(ivel,1,gp) \lambda
+shdy(ivel,3,gp) &
+shdy(ivel,1,gp) \lambda shdx(ivel,4,gp) & shdy(ivel,1,gp) \lambda
+shdy(ivel,4,gp) \\
+
+shdx(ivel,2,gp) \lambda shdx(ivel,3,gp) & shdx(ivel,2,gp) \lambda
+shdy(ivel,3,gp) &
+shdx(ivel,2,gp) \lambda shdx(ivel,4,gp) & shdx(ivel,2,gp) \lambda
+shdy(ivel,4,gp) \\
+
+shdy(ivel,2,gp) \lambda shdx(ivel,3,gp) & shdy(ivel,2,gp) \lambda
+shdy(ivel,3,gp) &
+shdy(ivel,2,gp) \lambda shdx(ivel,4,gp) & shdy(ivel,2,gp) \lambda
+shdy(ivel,4,gp) \\
+
+shdx(ivel,3,gp) \lambda shdx(ivel,3,gp) & shdx(ivel,3,gp) \lambda
+shdy(ivel,3,gp) &
+shdx(ivel,3,gp) \lambda shdx(ivel,4,gp) & shdx(ivel,3,gp) \lambda
+shdy(ivel,4,gp) \\
+
+\null & shdy(ivel,3,gp) \lambda
+shdy(ivel,3,gp) &
+shdy(ivel,3,gp) \lambda shdx(ivel,4,gp) & shdy(ivel,3,gp) \lambda
+shdy(ivel,4,gp) \\
+
+\null & \null &
+shdx(ivel,4,gp) \lambda shdx(ivel,4,gp) & shdx(ivel,4,gp) \lambda
+shdy(ivel,4,gp) \\
+
+\null & \null &
+\null & shdy(ivel,4,gp) \lambda shdy(ivel,4,gp) \\
+\end{array} \right ].
+\end{equation} }
+
+Now putting all this together you have enough information to finish your code project.
+
+\end{document}
\ No newline at end of file
More information about the cig-commits
mailing list