### (Deleted) Instability of nonlinear elastodynamics

68
views
0
4 weeks ago by
I played around with the undocumented elastodynamics demo and had some numerical issues. I changed the definition of the stress (for the hyperelastic Saint-Venant Kirchhoff model) but kept almost everything else the same. With this material model, the equations are nonlinear, so instead of forms 'a' and 'L', now I have 'F' and 'J' with 'u' being a Coefficient. Here are the changed parts in the ufl:
# Stress tensor
def sigma(r):
I = Identity(r.geometric_dimension())
F0 = I + grad(r)
E = 0.5*(F0.T*F0-I)
return lmbda*tr(E)*I + 2.0*mu*E

# Linear form
E1 = factor_m1*inner(r, u)*dx  \
+factor_d1*inner(r, u)*dx \
+(1.0-alpha_f)*inner(grad(r), sigma(u))*dx
E2 = factor_m1*inner(r, u0)*dx + factor_m2*inner(r, v0)*dx + factor_m3*inner(r, a0)*dx  \
+ factor_d1*inner(r, u0)*dx + factor_d2*inner(r, v0)*dx + factor_d3*inner(r, a0)*dx  \
- alpha_f*inner(grad(r), sigma(u0))*dx \
+ inner(r, f)*dx + (1.0-alpha_f)*inner(r, p)*ds(3) + alpha_f*inner(r, p0)*ds(3)
F = E1 - E2

# Jacobian
du = TrialFunction(element)
J  = derivative(F, u, du)
​

In the cpp part, the appropriate values are assigned to the forms 'F' and 'J'. Following the Cahn-Hilliard demo, I use the Newton method in each timestep to solve the problem. These are the modified lines of the cpp code:
  // Create solver Parameters
Parameters params("nonlinear_variational_solver");
Parameters newton_params("newton_solver");
newton_params.add("relative_tolerance", 1e-10);
params.add(newton_params);

// Start time stepping
std::size_t step = 0;
while (t < T)
{
// Update for next time step
t += *dt;
cout << "Time: " << t << endl;

// Solve nonlinear variational problem
solve(F == 0, *u, bc, J, params);

// Update velocity and acceleration
update_a(a, *u, *a0, *v0, *u0, *beta, *dt);
update_v(v, a, *a0, *v0, *gamma, *dt);
*u0 = *u; *v0 = v; *a0 = a;

Save solutions to file
if (step % 1 0 == 0)
{
file_u << *u;
}
++step;
}
​

The results seem to be accurate for some timesteps but after some time the Newton iteration fails. It seems to me, that the problem occurs when the deformation changes direction. I tried to reduce the timestep but the problem remained. What can be the reason for this instability? As far as I am concerned, the generalized-$\alpha$α method should work for nonlinear problems as well.

Here is the deformation of a cantilever beam demonstrating the problem:

I post the Newton iteration before the error and the error:
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 3.049e-02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-10)
Newton iteration 1: r (abs) = 6.863e-03 (tol = 1.000e-10) r (rel) = 2.251e-01 (tol = 1.000e-10)
Newton iteration 2: r (abs) = 5.976e-04 (tol = 1.000e-10) r (rel) = 1.960e-02 (tol = 1.000e-10)
Newton iteration 3: r (abs) = 1.934e-05 (tol = 1.000e-10) r (rel) = 6.342e-04 (tol = 1.000e-10)
Newton iteration 4: r (abs) = 4.702e-08 (tol = 1.000e-10) r (rel) = 1.542e-06 (tol = 1.000e-10)
Newton iteration 5: r (abs) = 3.160e-13 (tol = 1.000e-10) r (rel) = 1.036e-11 (tol = 1.000e-10)
Newton solver finished in 5 iterations and 5 linear solver iterations.
Time: 6.825
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 3.174e-02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-10)
Newton iteration 1: r (abs) = 8.671e-03 (tol = 1.000e-10) r (rel) = 2.732e-01 (tol = 1.000e-10)
Newton iteration 2: r (abs) = 1.097e-03 (tol = 1.000e-10) r (rel) = 3.456e-02 (tol = 1.000e-10)
Newton iteration 3: r (abs) = 2.184e-04 (tol = 1.000e-10) r (rel) = 6.883e-03 (tol = 1.000e-10)
Newton iteration 4: r (abs) = 9.723e-03 (tol = 1.000e-10) r (rel) = 3.064e-01 (tol = 1.000e-10)
Newton iteration 5: r (abs) = 2.484e-03 (tol = 1.000e-10) r (rel) = 7.826e-02 (tol = 1.000e-10)
Newton iteration 6: r (abs) = 6.795e-04 (tol = 1.000e-10) r (rel) = 2.141e-02 (tol = 1.000e-10)
Newton iteration 7: r (abs) = 2.479e-04 (tol = 1.000e-10) r (rel) = 7.809e-03 (tol = 1.000e-10)
Newton iteration 8: r (abs) = 4.503e-04 (tol = 1.000e-10) r (rel) = 1.419e-02 (tol = 1.000e-10)
Newton iteration 9: r (abs) = 2.143e-04 (tol = 1.000e-10) r (rel) = 6.752e-03 (tol = 1.000e-10)
Newton iteration 10: r (abs) = 2.080e-02 (tol = 1.000e-10) r (rel) = 6.552e-01 (tol = 1.000e-10)
Newton iteration 11: r (abs) = 5.253e-03 (tol = 1.000e-10) r (rel) = 1.655e-01 (tol = 1.000e-10)
Newton iteration 12: r (abs) = 1.369e-03 (tol = 1.000e-10) r (rel) = 4.313e-02 (tol = 1.000e-10)
Newton iteration 13: r (abs) = 4.055e-04 (tol = 1.000e-10) r (rel) = 1.278e-02 (tol = 1.000e-10)
Newton iteration 14: r (abs) = 2.144e-04 (tol = 1.000e-10) r (rel) = 6.754e-03 (tol = 1.000e-10)
Newton iteration 15: r (abs) = 1.833e-02 (tol = 1.000e-10) r (rel) = 5.775e-01 (tol = 1.000e-10)
Newton iteration 16: r (abs) = 4.636e-03 (tol = 1.000e-10) r (rel) = 1.461e-01 (tol = 1.000e-10)
Newton iteration 17: r (abs) = 1.215e-03 (tol = 1.000e-10) r (rel) = 3.828e-02 (tol = 1.000e-10)
Newton iteration 18: r (abs) = 3.686e-04 (tol = 1.000e-10) r (rel) = 1.161e-02 (tol = 1.000e-10)
Newton iteration 19: r (abs) = 2.193e-04 (tol = 1.000e-10) r (rel) = 6.911e-03 (tol = 1.000e-10)
Newton iteration 20: r (abs) = 2.149e-03 (tol = 1.000e-10) r (rel) = 6.770e-02 (tol = 1.000e-10)
Newton iteration 21: r (abs) = 5.965e-04 (tol = 1.000e-10) r (rel) = 1.879e-02 (tol = 1.000e-10)
Newton iteration 22: r (abs) = 2.324e-04 (tol = 1.000e-10) r (rel) = 7.322e-03 (tol = 1.000e-10)
Newton iteration 23: r (abs) = 7.238e-04 (tol = 1.000e-10) r (rel) = 2.281e-02 (tol = 1.000e-10)
Newton iteration 24: r (abs) = 2.568e-04 (tol = 1.000e-10) r (rel) = 8.091e-03 (tol = 1.000e-10)
Newton iteration 25: r (abs) = 3.830e-04 (tol = 1.000e-10) r (rel) = 1.207e-02 (tol = 1.000e-10)
Newton iteration 26: r (abs) = 2.167e-04 (tol = 1.000e-10) r (rel) = 6.827e-03 (tol = 1.000e-10)
Newton iteration 27: r (abs) = 4.022e-03 (tol = 1.000e-10) r (rel) = 1.267e-01 (tol = 1.000e-10)
Newton iteration 28: r (abs) = 1.062e-03 (tol = 1.000e-10) r (rel) = 3.346e-02 (tol = 1.000e-10)
Newton iteration 29: r (abs) = 3.324e-04 (tol = 1.000e-10) r (rel) = 1.047e-02 (tol = 1.000e-10)
Newton iteration 30: r (abs) = 2.328e-04 (tol = 1.000e-10) r (rel) = 7.335e-03 (tol = 1.000e-10)
Newton iteration 31: r (abs) = 7.113e-04 (tol = 1.000e-10) r (rel) = 2.241e-02 (tol = 1.000e-10)
Newton iteration 32: r (abs) = 2.542e-04 (tol = 1.000e-10) r (rel) = 8.010e-03 (tol = 1.000e-10)
Newton iteration 33: r (abs) = 3.992e-04 (tol = 1.000e-10) r (rel) = 1.258e-02 (tol = 1.000e-10)
Newton iteration 34: r (abs) = 2.148e-04 (tol = 1.000e-10) r (rel) = 6.769e-03 (tol = 1.000e-10)
Newton iteration 35: r (abs) = 1.068e-02 (tol = 1.000e-10) r (rel) = 3.367e-01 (tol = 1.000e-10)
Newton iteration 36: r (abs) = 2.726e-03 (tol = 1.000e-10) r (rel) = 8.588e-02 (tol = 1.000e-10)
Newton iteration 37: r (abs) = 7.394e-04 (tol = 1.000e-10) r (rel) = 2.330e-02 (tol = 1.000e-10)
Newton iteration 38: r (abs) = 2.600e-04 (tol = 1.000e-10) r (rel) = 8.193e-03 (tol = 1.000e-10)
Newton iteration 39: r (abs) = 3.653e-04 (tol = 1.000e-10) r (rel) = 1.151e-02 (tol = 1.000e-10)
Newton iteration 40: r (abs) = 2.201e-04 (tol = 1.000e-10) r (rel) = 6.936e-03 (tol = 1.000e-10)
Newton iteration 41: r (abs) = 1.896e-03 (tol = 1.000e-10) r (rel) = 5.975e-02 (tol = 1.000e-10)
Newton iteration 42: r (abs) = 5.343e-04 (tol = 1.000e-10) r (rel) = 1.683e-02 (tol = 1.000e-10)
Newton iteration 43: r (abs) = 2.226e-04 (tol = 1.000e-10) r (rel) = 7.015e-03 (tol = 1.000e-10)
Newton iteration 44: r (abs) = 1.393e-03 (tol = 1.000e-10) r (rel) = 4.390e-02 (tol = 1.000e-10)
Newton iteration 45: r (abs) = 4.114e-04 (tol = 1.000e-10) r (rel) = 1.296e-02 (tol = 1.000e-10)
Newton iteration 46: r (abs) = 2.141e-04 (tol = 1.000e-10) r (rel) = 6.745e-03 (tol = 1.000e-10)
Newton iteration 47: r (abs) = 3.521e-02 (tol = 1.000e-10) r (rel) = 1.109e+00 (tol = 1.000e-10)
Newton iteration 48: r (abs) = 8.856e-03 (tol = 1.000e-10) r (rel) = 2.791e-01 (tol = 1.000e-10)
Newton iteration 49: r (abs) = 2.269e-03 (tol = 1.000e-10) r (rel) = 7.149e-02 (tol = 1.000e-10)
Newton iteration 50: r (abs) = 6.262e-04 (tol = 1.000e-10) r (rel) = 1.973e-02 (tol = 1.000e-10)
terminate called after throwing an instance of 'std::runtime_error'
what():

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2017.2.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
Community: FEniCS Project

### 1 Answer

0
29 days ago by
Try using mumps. Also post complete error that your are facing. Are you using default solver ?
I use the default solver. Trying to set mums as a linear solver did not help unfortunately. I added the error into my question.
written 29 days ago by fehere
Please login to add an answer/comment or follow this question.
The thread is closed. No new answer/comment may be added.

Similar posts:
Search »
• Nothing matches yet.