9

I'm working in Mathematica to plot a 2-d catenary in Euclidean space according to this Math Stack Exchange Link, specified with a square Dirichlet boundary condition.

Here's what I've coded using Mathematica, followed by an error I receive. My goal was to set the surface integral to 1.5 because it's gotta be at least 1 if it's a sagging surface pinned at edges on a $1\times 1$ square.

NDSolve[
  { (-(D[u[x,y],x,x]+D[u[x,y],y,y])*(u[x,y]+λ)) == 1
  , Integrate[
      ( 1+D[u[x,y],x]+D[u[x,y],y] )^(1/2) ∈ RectangularPolygon[4],x,y]==1.5
    , DirichletCondition[u[x, y] == 0, True]
  }
, u
, {x, y} ∈ RegularPolygon[4]
]

NDSolve: There are fewer dependent variables, {u[x,y]}, than equations, so the system is overdetermined

I think this is conditionally stable, so I've worked to provide a fair approximation of a stable error. If I omit the surface integral line, I'm told that NDSolve doesn't deal with non-linear coefficients. I was told elsewhere that this nonlinearity error is a hint that more detail is needed.

OK, I took a stab at improving the translation. I've got a coefficient array error that NDSolve can't parse, but it's closer to the problem statement. Surface integral not included yet.

NDSolve[
{
Sqrt[1+Norm[Grad[u[x,y],{x,y}]]^2 ]*Div[Grad[u[x,y],{x,y}]/Sqrt[1+Norm[Grad[u[x,y],{x,y}]]^2 ] ,{x,y}]==1/(u[x,y]+1),
DirichletCondition[u[x, y] == 0, True]
}, 
u, 
{x, y} ∈ RegularPolygon[4]
]

image of chamois cloth

The trouble with this problem is that a physical analog is basically impossible. Any two-dimensional approximation of a heavy surface, such as a sheet of chain male, introduces anisotropy and also costs a lot. A sheet of chamois was used to get a first stab, but, as shown, it seems to have something different going on in the corners.

Work done by @xzczd below shows a shape more like the Laplacian over a similar domain, possessing the similar, albeit less-pronounced, re-entrant geometry along a geodesic from the peak to the vertex (see here; Google image result).

What I'll probably need to do is generate a .CSV file with this and then hand-prototype a form from which I can cast a concrete part. CNC foam is a nice option for a daydream, but it costs about \$1500 x 2 where I'm at in the US Midwest. Rather than shell out the cash on something that might melt at cone 10 (pottery term), I may be better off measuring a grid of sticks, draping plaster gauze cloth over the peaks and smoothing subsequent dollops of hydrocal thereover with an eye on a printout to capture detail. I feign no skill at anticipating a successful method here, as I've not worked on something this big before.

It could come in handy that this surface looks like the native "dome" feature offered in SolidWorks, but, without sufficient documentation regarding the algorithm therein used, and also with there already being discrepancies between the Laplacian already looking 10% off from the @xzczd, I'm likely to err on the side of extra elbow grease and precision. If we've got the right math here, why use anything else?

A note on the picture above; chamois cloth (used wet as shown) stretches as it's comprised of, basically, an epidermal sponge, so there's probably some heat energy (and a smidge of elastic energy) absorbing the gravitational potential that deprives the final shape of any proximity to a more truthful catenary surface.

Dope work, all. Thanks to @xzczd.

Foam board

@xzczd - This is a 38% scale prototype I'm making. If I can find affordable CNC service I'll use it, otherwise I'll make a pantograph and scale this up using a drill holding a ball mill following a ball stylus that runs over this form like a cam/template.

The ribs shown were generated using your Mathematica code modified a bit to output a $10 \times 10$ grid of points. I'm working with half of the shape at first (which is what you see). The next step will be to add the remaining webbing-like shapes that I've made from card stock.

Thank you - I think this is the first time in the history I've learned that I've seen this shape be made. I hope it works.

Been seeing similar shapes on domes for a few weeks now so doubt the immediate above.

Progress pics below. One shows a gauze cloth attempt that failed, but the next shows how I got to make the slip cast mold from solid plaster (similar in concept to the gauze cloth but without the cloth—two thin shells were made by pouring nearly-cured plaster slurry over the yellow papier mache and polyurethane master, then those two shells were joined with another thick overlay of plaster).

Current state, leather hard. Slip cast the main body and added hand built clay handles. Drying slowly now then firing and then testing. God help us all with the picture upload.

It’s probably 50-80% accurate. One big issue will be that the analytical surface is what I started with for the mold—the final shell thickness isn’t built in. If I get another chance I’ll do it, but for this moment I just want to know if I can cast something this large and then transport it to the field in one piece. If not, I may have to use castable refractory or try another approach.

enter image description here

enter image description here

enter image description hereenter image description here

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Brandon
  • 93
  • 6
  • 1
    Your translation for the PDE in the link is wrong, the $\sqrt{1+|\nabla u|^2}$ term cannot be cancelled, notice one of them is inside $\nabla\cdot{}$. – xzczd Jun 21 '18 at 11:26
  • @xzczd: I'll work to transcribe the PDE as-written to make this more fun to follow in the future. – Brandon Jun 21 '18 at 15:11
  • @Mariusz: Lambda is a Lagrange multiplier; rectangular polygon herein specified with 4 vertices is a cute way to require a square domain (similar to Disk). – Brandon Jun 21 '18 at 15:15
  • Do you mean that, you just want to solve this equation as a simplified one at the moment? – xzczd Jun 21 '18 at 15:32
  • @mariusz probably means there's no built-in function named as RectangularPolygon… do you mean RegularPolygon at that position? – xzczd Jun 21 '18 at 17:17
  • @xzczd I did mean it and took a stab at it - see edited post; looks like more folks have chimed in with more skill, but at least I learned how to properly use Div and Grad. – Brandon Jun 21 '18 at 20:40
  • @mariusz - RegularPolygon works in other problems - documentation here http://reference.wolfram.com/language/ref/RegularPolygon.html – Brandon Jun 21 '18 at 21:52
  • The 2D case is more sophisticated as the 1D case because one cannot work with isometric deformations. I don't think that the model elaborated in https://mathoverflow.net/a/69837 is physically meaningful. In particular, the area constraint is nonsense: It is precisely the (local!) change of area that produces the stretch forces that counter the gravitational forces. So it's no wonder that the physical experiments do not fill well to simulated data. – Henrik Schumacher Aug 15 '18 at 07:30
  • IMHO, one has to add a stretching term to the total energy of the system and the gravitational energy term should read $\int_\varOmega u , \operatorname{d} !x$ (it has to involve only the area element of the surface in rest state; it must not dependend of the area element of the deformed surface). – Henrik Schumacher Aug 15 '18 at 07:32
  • @Henrik, I couldn’t come up with a 2D analog of what a 1D chain does. A gold necklace is slinky, and the changes in distance between hinge fulcrums is generally fixed as compared to those of a sheet of chain mail, local portions of which will go from square to trapezoidal under shear. No shear should occur in the weight loading of a 2D catenary some (i.e., by magically bringing it from outer space to the surface of a planet flatly). The chamois cloth then stretches, which a necklace (or other chain) doesn’t do when it’s used to trace a catenary for a masonry arch. – Brandon Aug 16 '18 at 14:17
  • @xzczd uploaded pic – Brandon Aug 21 '18 at 17:53

1 Answers1

8

I won't say it's easy to solve.

"I'm told that NDSolve doesn't deal with non-linear coefficients. " Yes, currently NDSolve cannot deal with nonlinear stationary PDE, so let's turn to finite difference method (FDM) to discretize $$\sqrt{1+|\nabla u|^2}\ \nabla\cdot{}\frac{\nabla u}{\sqrt{1+|\nabla u|^2}}=\frac{1}{u+\lambda}$$. I'll use pdetoae for the generation of difference equations.

To take the constraint $$\int_{\Omega}\sqrt{1+|\nabla u|^2}dx=A$$ into consideration, I'll use trapezoidal rule to discretize it to a algebraic equation and solve it together with those difference equations.

points = 24;
{bL, bR} = domain = {0, 1};
grid = Array[# &, points, domain];
grad = Grad[#, {x, y}] &;
div = Div[#, {x, y}] &;
A = 3/2;
With[{u = u[x, y]}, sqrt = Sqrt[1 + grad[u].grad[u]]; 
  eq = sqrt div[grad[u]/sqrt] (u + λ) == 1;
  {bc@x, bc@y} = Table[u == 0 /. var -> b, {var, {x, y}}, {b, domain}]];
sqrtfunc = Function[{x, y}, #] &@sqrt;
difforder = 2;
(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
ptoafunc = pdetoae[u[x, y], {grid, grid}, difforder];
del = #[[2 ;; -2]] &;
ae = del /@ del@ptoafunc[eq];
aebc@x = del /@ ptoafunc[bc@x];
aebc@y = ptoafunc[bc@y];
rule = Flatten /@ (Outer[sqrtfunc, grid, grid] -> ptoafunc[sqrt]) // Thread;
trap[value_, grid_: grid] := 1/2 Differences[grid] (Most@value + Rest@value) // Total;
int = trap[Function[y, trap[Function[x, sqrtfunc[x, y]] /@ grid]] /@ grid] /. rule;
solfuncguess[x_, y_] := 0; λguess = 1/3;
solrule = FindRoot[{ae, aebc /@ {x, y}, int == A}, 
    Append[Flatten[#, 1] &@
      Table[{u[x, y], solfuncguess[x, y]}, {x, grid}, {y, 
        grid}], {λ, λguess}], MaxIterations -> 1000]; // AbsoluteTiming
(* {35.472657, Null} *)
solfunc = ListInterpolation[
   Partition[Most[solrule][[All, -1]], points], {domain, domain}];
solrule // Last
(* λ -> 0.555491 *)

NIntegrate[sqrt /. u -> solfunc, {x, bL, bR}, {y, bL, bR}]
(* 1.51602 *)
Plot3D[{solfunc[x, y]}, {x, 0, 1}, {y, 0, 1}]

Mathematica graphics

Remark

  1. The achieved $A$ value isn't perfect, to achieve a better $A$, better initial guess is probably needed, which seems not to be that easy to find for this problem.

  2. Guys read my relevant posts before (Well, will there be any?) may feel surprised that I'm not using Gauss-Legendre quadrature formula for numeric integration this time, that's because, with Gauss-Legendre quadrature formula, the solution never converges to a reasonable position. I don't know why.

  3. Honestly speaking I'm not 100% sure if the solution above is reliable.

Community
  • 1
xzczd
  • 65,995
  • 9
  • 163
  • 468