Numerical modelling and simulations of nonlinear elasto-plasticity for geodynamics applications


Monday, 12 December, 2016 - 17:30

Speaker: Mr Roberto Porcù, Politecnico di Milano, Department of Mathematics, MOX

Room: A-005, SISSA Campus, Via Bonomea 265, Trieste

The aim of this work is the numerical simulation of the Earth's lithosphere extension, with a particular focus on continental rifting that is a tectonical process which has a deep influence on other geodynamical phenomena and on the spatial distribution of some important natural resources.  During last decades a large amount of research concerning this topic has been performed by means of analogical and numerical models, named sandboxes.  This work is motivated by the observation of several issues inherent to the analogic experiments and by the generally accepted adoption, in the known literature, of purely linear rheologies for numerical tests, which are clearly incoherent with the natural finite-strain deformation regime of the problem at hand.  In this work a fully-nonlinear elasto-plastic rheological approach has been developed and implemented. The variables of the problem are parametrized by material coordinates thus a totally-Lagrangian formulation has been adopted together with a quasi-static motion assumption justified by the analysis of the dimensionless quantities of the governing equations.  In order to avoid the volumetric locking arising from the local strong enforcement of a nearly incompressible deformation, a three-fields weak variational formulation of the problem has been taken into account, as discussed in [Simo1998]: this method consists in the weak enforcement of the volumetric constraint by introduction of a new independent dilatation variable together with its dual variable that can be interpreted as the Kirchhoff pressure.  This leads to an additive decomposition of the Helmoltz free energy into its isochoric and volumetric parts.  Coherently to the fully-nonlinear rheological approach the deformation gradient is multiplicatively decomposed into its elastic and plastic parts, according to the Kroner-Lee decomposition.  The global nonlinear system is solved by application
of the Newton-Raphson linearization method. Then the constitutive update, which is computed point-wise on the quadrature nodes, is carried out firstly determining an elastic trial state; if the trial state is nonadmissible, the exponential return map algorithm is performed in order to recover admissibility. The nonlinear return map scheme is solved by means of the Newton-Raphson linearization method.  A Drucker-Prager yield criterion and a non-associative flow rule have been considered since they represent really good choices for geological applications, as stated in [BelytschkoLiu2001].  A new C++
MPI-parallel code has been implemented; it is based on the deal.II finite element library; it handles three-dimensional multi-layered domains.  The results are encouraging since they show the natural development of V-shaped shear bands which is a deformation pattern that the most-diffused softwares in the geodynamics field are able to catch only if a weak seed is artificially added into the domain.