Problem Set 4 Group

Problem G4-1

It turns out that when the materials properties grouping = 1, this simplifies future calculations tremendously and allows the graphs to look a little prettier.  However, a realistic Poisson's ratio has been chosen, as well as arbitrary G, B.

G  =  1;
b  =  1;
ν  =   .3;
MaterialProp  =  
        G b/(2 Pi (1 - ν));
\(\(StressMatrixXY\  = \ 
        MaterialProp\ {
				{\((\(\(-y\)\ \((3  x\^2 + y\^2)\)\)\/\((x\^2 + y\^2)\)\^2)\), \ \((\(x\ \((x\^2 - y\^2)\)\)\/\((x\^2 + y\^2)\)\^2)\)}, 
				{\ \((\(x\ \((x\^2 - y\^2)\)\)\/\((x\^2 + y\^2)\)\^2)\), \((\(y\ \((x\^2 - y\^2)\)\)\/\((x\^2 + y\^2)\)\^2)\)}
				};\)\)
StressStateXY  = StressMatrixXY;
StressStateXY // MatrixForm

a.

PolarConditions  =  { x -> r Cos[θ],  
          y -> r Sin[θ]};
StressStateRθ  =  
        StressStateXY /. PolarConditions;
StressStateRθ // Simplify // MatrixForm

b.

The hydrostatic pressure should be an invariant under any axial manipulation, and can be calculated by summing the diagonal.

HydroPressXY  =  (1/3) (StressStateXY[[1, 1]] + 
              StressStateXY[[2, 2]]);
HydroPressRθ  = (1/3) ( 
            StressStateRθ[[1, 1]] + 
              StressStateRθ[[2, 
                  2]]);
HydroPressRθ // Simplify
 HydroPressXY // Simplify

Plot3D[HydroPressXY,  {x, -1, 1}, {y, -1, 1},  
        ViewPoint -> {1, 1, 1}];

And now a plot of the isobars.

ContourPlot[HydroPressXY,  {x, -1, 1}, {y, -1, 1}]

c.

Note for Part (c).  If you replace " := " with " = ",  you are inserting immediate assignments for the following delayed assignments.  This unfortunately can render two plots of the same equations very different.  The solution using immediate assignment (using "=") produces many errors when calculating the angle rotation because of denominators approaching 0, so I've chosen to present a solution utilizing the delayed assignment " := " instead.  What's most important in the solution, however, is understanding the process involved in solving for the angle and shear stress in the next two parts.

StressStateXY

UnitEigenVec := Eigenvectors[StressStateXY][[1]]/
        Norm[Eigenvectors[StressStateXY][[1]]]

CosofAngle := UnitEigenVec . {1, 0}

Plot3D[ArcCos[CosofAngle], {x, -10, 10}, {y, -10, 10}, 
        PlotPoints -> 100];

d.

Maximum shear stress vs. position...  From the Mohr circle construction, we know that the maximum shear stress at a point is just the radius of the Mohr circle, and this occurs at a 2θ angle of π/2 (i.e., a θ angle of π/4) with respect to the principal axes.
The radius of the Mohr circle is equal to (σmax-σmin)/2 (half the difference of the principal normal stresses), and σmax and σmin are simply obtained as the eigenvalues of the stress matrix.

MaxShearStress  :=  (Eigenvalues[
                StressStateXY][[1]] - Eigenvalues[
                StressStateXY][[2]])/2

Plot3D[MaxShearStress,  {x, -10, 10}, {y, -10, 10},  
      PlotPoints  -> 100]

Problem G4-2

a.

Parametric representation of a catenoid.  Though it was discussed in recitation that parametric equations for the catenary would suffice, here are the parametric equations for the 3D catenoid.  Typically, the variables u & v are used as parametric variables in more than one dimension (t is used for one dimensional parameterization).

X1  =  A Cosh[(v + B)/A] Sin[u]
Y1  =  A Cosh[(v + B)/A] Cos[u]
Z1 =  v

The catenary is simply:

X2 = t
Y2  =  A Cosh[(t + B)/A]

b.

The curvature of the catenary can be calculated as follows:

Clear[X, Y, A, B]

Y  =   A Cosh[(X + B )/A]

The radius of curvature of a function is defined as ((\((1 + y'\^\(\(\ \ \)\(2\)\))\))^(3/2))/ y''.  The "curvature" is simply the inverse of the radius of curvature.

RadCurvature  =  (((1  +  (((D[Y, X]))^2)))^(3/2))/(D[
              D[Y, X], X]);
Curvature  =  1/RadCurvature

First a plot of the original function.

Plot[Y, {X, -5, 5}]

And now a plot of the curvature.

Plot[Curvature,  {X, -5, 5}]

c.

This part is pretty self-explanatory.  Plug into Mathematica and plot.

PartC  =  Curvature  -  1/Y;

Plot[PartC,  {X, -5, 5}];

d.

Note, the following two parts were edited from the solution turned in by the Asphodel group, consisting of Yun, Handler, Cooper, Gray, Duffy.  Their approach was methodical and logical.

\(Take\ R\_1 = 
      R\_2 . \ \ One\ would\ presume\ that\ R\ and\ L\ scale\ the\ catenary . \ \ A\ natural\ parameter\ is\ \(L\/R\) . \ \ Parameters\ are\ usully\ \(\(dimensionless\)\(.\)\)\)

Catenary  = 
      y - A Cosh[((x + B)/A)] /. {y -> (yp  R), 
          x -> (xp R), A -> (ap  R), 
          B ->  (bp R)}

\(Taking\ y \((0)\)\  = \ \(R\ and\ y \((L)\) = R\), \ 
    we\ get\ simultaneous\ equations\ in\ a\ and\ b\ that\ can\ be\ solved\ for\ a\ given\ L\/R\)

Editor's Note: The following several steps are the nuts and bolts of the above described algebra.

Catenary /. {yp -> R, xp -> 0}

FullSimplify[((R^2) - ap R Cosh[(bp/ap)])/R]

FullSimplify[(Catenary /. {yp -> R, xp -> L})/
        R]

Catenary = 
      y - A Cosh[((x + B)/A)] /. {A -> (ap R), 
          B ->  (bp R)}

Equation1 = Catenary /. {y -> R, x -> 0}

FinalEquation1 = FullSimplify[Equation1/R]

Equation2 = Catenary /. {y -> R, x -> L} /. 
        L -> R 1.2

FinalEquation2 = 
      1 - ap Cosh[(L/(ap R)) + ((bp R)/(ap R))]

Editor's Note: A function was defined that took an L and R value and returned values for A and B, as well as the catenary itself.

Given[L_, R_] := 
      Plot[ap R Cosh[((bp R + x)/(ap R))] /. 
            FindRoot[{0 == 1 - ap Cosh[(bp/ap)], 
                0 == 1 - ap Cosh[(bp/ap) + (L/(ap R))]}, {ap,  .3, 
                5}, {bp,  .5, 8}, MaxIterations -> 1000], {x, 0, L}, 
          PlotRange -> {0, R}] && 
        Print["A is ", 
          ap R /. FindRoot[{0 == 1 - ap Cosh[(bp/ap)], 
                0 == 1 - ap Cosh[(bp/ap) + (L/(ap R))]}, {ap,  .3, 
                5}, {bp,  .5, 8}], " and B is ", 
          bp R /. FindRoot[{0 == 1 - ap Cosh[(bp/ap)], 
                0 == 1 - ap Cosh[(bp/ap) + (L/(ap R))]}, {ap,  .3, 
                5}, {bp,  .5, 8}]]

Check out this cool stuff!  It is clear that as you pull the rings farther apart, there is a larger dip in the soap bubble.  Note that for some ratios of L to R, Mathematica has difficulty converging while finding roots, which we will explain below.

Given[12, 18]

Given[20, 18]

Editor's Note: Note the appropriate boundary conditions.

We know that y/A=Cosh[(B+x)/A].  Let's take the case where x=0:R/A=Cosh[B/A].  We know that (ap R)=A and (bp R)=B.  Substituting this in to the previous equation, we get 1/ap=Cosh[bp/ap] Because bp and ap are tied to L and R,the ratio of bp to ap determines if solutions exist and how many solutions exist.  We can look at the graphs and determine how many times they intersect with ap as the variable.

SecondGraph = 
      Plot3D[Cosh[bp/(ap)], {ap, 0, 2}, {bp, 0, 1}, 
        Mesh -> False]

FirstGraph = 
      Plot3D[1/(ap), {ap, 0, 2}, {bp, 0, 1}, 
        Mesh -> False]

One can see the intersection of the graphs.  Where this occurs, there is an allowable ratio of L to R. Where there is no intersection, the soap bubble does not exist.

Show[FirstGraph, SecondGraph]

Plotting 1/ap and Cosh[bp/ap] on the same 2-D graph, we can see the allowable ratio bp/ap for different bp's- when there are 2 solutions, 1 solution, No solutions.

Plot[{Cosh[2/(ap)], 1/ap}, {ap, 0, 2}, 
      PlotRange -> {0, 5}]

Plot[{Cosh[1/(ap)], 1/ap}, {ap, 0, 2}, 
      PlotRange -> {0, 5}]

The ratio of 2/3 is very close to having only one solution

Plot[{Cosh[ .67/(ap)], 1/ap}, {ap, 0, 1.5}, 
      PlotRange -> {0, 6}]

Plot[{Cosh[ .5/(ap)], 1/ap}, {ap, 0, 1.5}, 
      PlotRange -> {0, 6}]

Plot[{Cosh[ .25/(ap)], 1/ap}, {ap, 0, 1.5}, 
      PlotRange -> {0, 20}]

e.

Here's a function that finds the equation of a catenary for a given L and R:

CatCalc[L_, R_] := 
      ap R Cosh[((bp R + x)/(ap R))] /. 
        FindRoot[{0 == 1 - ap Cosh[(bp/ap)], 
            0 == 1 - ap Cosh[(bp/ap) + (L/(ap R))]}, {ap,  .3, 
            5}, {bp,  .5, 8}]

CatCalc[18, 19]

Here we define a function that defines the surface area for a given R and L:

Area[L_, R_] := 
      Integrate[
        2 Pi CatCalc[L, R] (Sqrt[1 + D[CatCalc[L, R], x]^2]), {x, 0, 
          L}]

Area[18, 19]

\(We\ took\ E \((R, 
        L)\)\ to\ be\ equal\ to\ the\ Area\ times\ the\ surface\ tension . \ \ Therefore, \ \(E \((R, L)\)\)\/γRL = Area\/RL\)

Plot3D[Area[L, R]/(L R), {L, 18, 20}, {R, 18, 20}, 
      PlotPoints -> 15]