# Nonlinear elasticity PDE in Mathematica 12

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.

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.

• (+1), it's usually better to solve time dependent problems with the method of lines and FEM/TGP as a spatial discretization then as a pure spatial problem like done here. Usually specifying Method->"MethofOfLInes" then indicates some setup problem for the time dependent problem and is ultimately the reason a pure spatial FEM as chosen. – user21 Dec 19 '19 at 9:50
• So, for a problem like this, would it be best to use NDSolve only for the spatial problem (say, with default properties, like in my ODE example) and then write a loop (say simple Euler scheme) iterating over time? Or is there a magic NDSolve trick that could specify different methods for spatial and temporal integration as options? – Alexander Erlich Dec 26 '19 at 14:59
• @AlexanderErlich "……Euler scheme……" I doubt, because we need to transform the system to the standard form required by Euler scheme, which is not easy for the specific problem. (Related: mathematica.stackexchange.com/a/158519/1871) “…magic NDSolve trick…" There doesn't seem to be, because NDSolve is having difficulty in noticing MethodOfLines can be used in t direction. (Adding Method -> {MethodOfLines, TemporalVariable -> t} causes tvic warning. ) – xzczd Dec 27 '19 at 4:15
• @AlexanderErlich Anyway, I've added a solution for the toy example forcing NDSolve to use DAE solver in t direction. It'll probably fail for more complicated problem, though. – xzczd Dec 27 '19 at 5:17