# Solving system of nonlinear PDEs

My aim is to solve a system of nonlinear PDEs arising in nonlinear elasticity. I am new to Mathematica so I started by modifying this example. I tried to change it to neo-Hookean solid. The resulting system I am trying to solve looks like this:

\begin{aligned} \mathrm{div} \text{ } \mathbb{T}_R &= \mathbb{0} \\ \mathbb{T}_R &= 2 C (\mathbb{F} - \mathbb{F}^{-T}) + 2 D (\det \mathbb{F} - 1) (\det \mathbb{F}) \mathbb{F}^{-T} \\ \mathbb{F} &= \nabla \begin{pmatrix} u(x, y) \\ v(x, y) \end{pmatrix} \end{aligned}

where $$u$$ and $$v$$ are unknown. I arrived to the following code:

def[x_, y_] = {u[x, y], v[x, y]};
F[x_, y_] = Inactive[Grad][def[x, y], {x, y}];

C = 1;
D = 1;

T[x_, y_] = 2*C*F[x, y] + 2*D*((Det[F[x, y]] - 1) Det[F[x, y]] - 2*C)*
Inverse[Transpose[F[x, y]]];

Subscript[Γ, D] = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x == 0];

op = Inactive[Div][T[x, y], {x, y}]

{ufun, vfun} = NDSolveValue[{op == {0, NeumannValue[-1, x == 5]},
Subscript[Γ , D]}, {u, v}, {x, 0, 5}, {y, 0, 1}];


But when I try to run this code, I get the following error messages:

NDSolveValue::underdet: There are more dependent variables, {u[x,y],v[x,y]}, than equations, so the system is underdetermined.
Set::shape: Lists {ufun,vfun} and NDSolveValue[…] are not the same shape.


• Don't use C and D as constant since they have another meanning in Mathematica. Apr 18, 2022 at 5:06
• The root problem, I believe, is that F appears nonlinearly in the PDE. In general, the finite element method as implemented in Mathematica, does not work well for highly nonlinear PDEs. The system also has a boundary condition problem at x = 0, where Inverse[Transpose[F[x, y]]] is singular. Apr 19, 2022 at 1:20
• @bbgodfrey I think any finite element solver has problems with this, this is nothing specific for the Mathematica FEM solver. Or can you suggest a solver that solves the problem given? Apr 19, 2022 at 5:06
• @user21 You are correct. Sorry for the poor choice of words. By the way, nice solutions below. (+1). Apr 19, 2022 at 13:14
• @bbgodfrey, no worries, I just wanted to check if I missed an important detail. Apr 19, 2022 at 13:44

There are some issue you'd need to look at, like material parameters and a singularity at 0 but here is how this would work in principal. You have a choice. Either you repetitively solve a stationary problem and approach the final value of the deformation or you convert this into a pseudo time dependent problem and let the time integrator approach the final configuration. In both cases some of your code needs to be fixed.

Static version:

def[x_, y_] = {u[x, y], v[x, y]};
F[x_, y_] = Grad[def[x, y], {x, y}];
C1 = 1;
D1 = 1;
detF = Det[F[x, y]];
T[x_, y_] =
2*C1*F[x, y] +
2*D1*((detF - 1) detF - 2*C1)*
Transpose[(detF*Inverse[F[x, y]])/(detF + $MachineEpsilon)]; op = Inactive[Div][#, {x, y}] & /@ T[x, y]; pfun = ParametricNDSolveValue[{-op == {0, NeumannValue[-p, x == 5]}, DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x == 0] }, {u, v}, {x, 0, 5}, {y, 0, 1}, p];  The detF +$MachineEpsilon is to work around the singularity. To show how this works in principal. Also note that you need a divergence of every equation. I also removed the Inactive[Grad] such that all coefficients are in the inactive div and will essentially be parsed as a DerivativePDETerm. Also check the sign for you equation of motion. You then start an iterative scheme with

pfun[0.000001]


And approach your final value of p.

Transient version: Here the idea is to let the time integrator do the work. Again, this shows how this works in principal, you'd need to set the problem up properly, material parameters etc.

def[x_, y_] = {u[t, x, y], v[t, x, y]};
F[x_, y_] = Grad[def[x, y], {x, y}];

C1 = 1;
D1 = 1;
rho = 1000;

detF = Det[F[x, y]];

T[x_, y_] =
2*C1*F[x, y] +
2*D1*((detF - 1) detF - 2*C1)*
Transpose[(detF*Inverse[F[x, y]])/(detF + \$MachineEpsilon)];

op = Inactive[Div][#, {x, y}] & /@ T[x, y];
Monitor[NDSolveValue[{rho*D[def[x, y], {t, 2}] - op == {0,
NeumannValue[-t*1, x == 5]},
DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, x == 0],
u[0, x, y] == 0, v[0, x, y] == 0,
Derivative[1, 0, 0][u][0, x, y] == 0,
Derivative[1, 0, 0][v][0, x, y] == 0
}, {u, v}, {t, 0, 1}, {x, 0, 5}, {y, 0, 1},
EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], monitor]


Additionally to what I said above I have added a density rho and made the initial conditions zero. This will start the time integration.

Keep in mind that solid mechanics problems are highly nonlinear and it is very unlikely that you will find a solution in one step - and that is not a problem with the Mathematica implementation but a general issue of nonlinear solid mechanics like hyperleastic material models with large deformations.