Nonlinear elasticity PDE in Mathematica 12

Question

Mathematica 12, Windows 10.

I am trying to solve a PDE in one spatial dimension $R$ and time $t$. I need a solution for displacement $r(R,t)$, radial Cauchy stress $T_{RR}(R,t)$, and radial growth $\gamma_{RR}(R,t)$.

To make things as simple as possible, I first write the static version of the system (two coupled ODEs), making sure we get a solution. Then I add a trivial time evolution turning it into a PDE, and hopefully we can figure out the problem together.

Static problem (coupled ODEs)

 γRR[R_, t_] := 3 
γθθ[R_, t_] := 4

eqns = {
  r[R, t] D[r[R, t], R] == γRR[R, t] γθθ[
     R, t] R,
  r[R, t]^4  R γθθ[R, t] D[TRR[R, t], R] == 
   2 γRR[R, 
     t] (r[R, t]^4 - R^4 γθθ[R, t]^4) 
  }

ic = {r[0.1, t] == 0.1, TRR[1, t] == 0}

sol = First@NDSolve[{eqns, ic}, {r[R, t], TRR[R, t]}, {R, 0.1, 1}]

PlotSol[sol_] := 
 Plot[{r[R, t], TRR[R, t]} /. sol /. t -> 1 // Evaluate, {R, 0.1, 1}, 
  AxesLabel -> {"R", None}, 
  PlotLegends -> Placed[{"r(R)", "TRR(R)"}, {Center, Top}]]
PlotSol@sol

The solution is as expected:

Dynamic problem (PDE)

Now I simply add a trivial equation to the system of differential equations: The time derivative of $\gamma_{RR}$ is zero, that is $\dot{\gamma}_{RR}=0$. I amend the initial conditions with $\gamma_{RR}(R, t=0)=3$, which of course completely decouples and should give the same result as in the ODE example. I now NDSolve now over a domain in both $R$ and $t$:

Clear@γRR

AppendTo[eqns, D[γRR[R, t], t] == 0]

AppendTo[ic, γRR[R, 0] == 3 ]

sol = First@
  NDSolve[{eqns, ic}, {r[R, t], TRR[R, t]}, {R, 0.1, 1}, {t, 0, 1}]
PlotSol@sol

However, unexpectedly, I get some error messages from the NDSolve and the correct solution is not returned:

What is going on here? I added a trivial uncoupled set of equations which should not have changed the results.

Answer 1

It's because NDSolve is using different method for these 2 problems. For the first problem, Shooting method, which is for nonlinear boundary value problem of ordinaray differential equation, is used; while for the second problem, FiniteElement method, which supports nonlinear problem as of v12, is used. Different solvers of NDSolve are not equally powerful. For the second problem, we can obtain a result by adjusting the option a little. (The result is not perfect though):

sol = First@
  NDSolve[{eqns, ic}, {r[R, t], TRR[R, t], γRR}, {R, 0.1, 1}, {t, 0, 1},
          InitialSeeding -> {r[R,t] == 3, TRR[R,t] == 3, γRR[R,t] == 2},
          Method -> {FiniteElement, MeshOptions -> MaxCellMeasure -> 0.05}]
PlotSol@sol

"Then, is it possible to avoid using pure FiniteElementin NDSolve?" Well, for this specific toy problem, if we differentiate the first 2 equations in t direction and add proper initial condition, it's possible:

γθθ[R_, t_] := 4
tend = 1;
bleft = 1/10; bright = 1;
With[{r = r[R, t], TRR = TRR[R, t], γRR = γRR[R, t], γθθ = γθθ[R, t]},

  eqns = {r D[r, R] == γRR γθθ R,
          r^4 R γθθ D[TRR, R] == 2 γRR (r^4 - R^4 γθθ^4),
          D[γRR, t] == 0};

  bc = {r == 0.1 /. R -> bleft, TRR == 0 /. R -> bright};

  ic = γRR == 3 /. t -> 0;];

Block[{γRR}, γRR[R_, t_] := 3; 
  solinitlst = 
   NDSolveValue[{eqns // Most, bc}, {r[R, t], TRR[R, t]}, {R, bleft, bright}]];

icadd = {r[R, 0], TRR[R, 0]} == solinitlst;

sollst = NDSolveValue[{D[eqns // Most, t], eqns // Last, bc, ic, icadd}, {r, 
    TRR, γRR}, {t, 0, tend}, {R, bleft, bright}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 50, 
       "MinPoints" -> 50, "DifferenceOrder" -> 2}}];

Plot[sollst[R, tend] // Through // Evaluate, {R, bleft, bright}]

Notice this approach will probably fail in more complicated cases, because the DAE solver NDSolve is not that strong. In those cases, advanced technique is needed, here is an example.