An upwind cell centred finite volume method for large strain explicit solid dynamics in OpenFOAM
Authors:   , Haider Jibran 
Published: July 2023

In practical engineering applications involving extremely complex geometries, meshing typically constitutes a large portion of the overall design and analysis time. In the computational mechanics community, the ability to perform calculations on tetrahedral meshes has become increasingly important. For these reasons, automated tetrahedral mesh generation by means of Delaunay and advancing front techniques have recently received increasing attention in a number of applications, namely: crash impact simulations, cardiovascular modelling, blast and fracture modelling. Unfortunately, modern industry codes in solid mechanics typically rely on the use of traditional displacement based Finite Element formulations which possess several distinct disadvantages, namely: (1) reduced order of convergence for strains and stresses in comparison with displacements; (2) high frequency noise in the vicinity of shocks; and (3) numerical instabilities associated with shear locking, volumetric locking and pressure checker-boarding.

In order to address the above mentioned shortcomings, a new mixed-based set of equations for solid dynamics formulated in a system of first order hyperbolic conservation laws was introduced. Crucially, the new set of conservation laws has a similar structure to that of the well known Euler equations in the context of Computational Fluid Dynamics (CFD). This enables us to borrow some of the available CFD technologies and to adapt the method in the context of solid dynamics. This thesis builds on the work carried out by Lee et al. 2013 by further developing the upwind cell centred finite volume framework for the numerical analysis of large strain explicit solid dynamics and its tailor-made implementation within the open source code OpenFOAM, extensively used in industrial and academic environments. The object oriented nature of OpenFOAM implementation provides a very efficient platform for future development. In this computational framework, the primary unknown variables are linear momentum and deformation gradient tensor of the system. Moreover, the formulation is further extended for an additional set of geometric strain measures comprising of the co-factor of deformation gradient tensor and the Jacobian of deformation, in order to simulate polyconvex constitutive models ensuring material stability.

The domain is spatially discretised using a standard Godunov-type cell centred framework where second order accuracy is achieved by employing a linear reconstruction procedure in conjunction with a slope limiter. This leads to discontinuities in variables at the cell interface which motivate the use of a Riemann solver by introducing an upwind bias into the evaluation of numerical contact fluxes. The acoustic Riemann solver presented is further developed by applying preconditioned dissipation to improve its performance in the near incompressibility regime and extending its range to contact applications. Moreover, two evolutionary frameworks are proposed in this study to satisfy the underlying involutions (or compatibility conditions) of the system. Additionally, the spatial discretisation is alternatively represented through a nodal cell centred finite volume framework for comparison purposes. From a temporal discretisation point of view, a two stage Total Variation Diminishing Runge-Kutta time integrator is employed to ensure second order accuracy. Additionally, inclusion of a global posteriori angular momentum projection procedure enables preservation of angular momenta of the system. Finally, benchmark numerical examples are simulated to demonstrate various aspects of the formulation including mesh convergence, momentum preservation and the locking-free nature of the formulation on complex computational domains.

En aplicaciones prácticas de ingeniería que implican geometrías extremadamente complejas, el mallado requiere típicamente una gran parte del tiempo total de diseño y análisis. En la comunidad de mecánica computacional, la capacidad de realizar cálculos sobre mallas tetraédricas está siendo cada vez más importante. Por estas razones, la generación automatizada de mallas tetraédricas por medio de técnicas de Delaunay y frente avanzado han recibido cada vez más atención en ciertas aplicaciones, a saber: simulaciones de impacto, modelado cardiovascular, modelado de explosión y fractura. Por desgracia, los códigos en la industria moderna para mecánica de sólidos se basan normalmente en el uso de formulaciones tradicionales de Elementos Finitos formulados en desplazamientos que poseen varias desventajas: (1) menor orden de convergencia para tensiones y deformaciones; (2) ruido de alta frecuencia cerca de las ondas de choque; y (3) inestabilidades numéricas asociadas con el bloqueo a cortante, el bloqueo volumétrico y oscilaciones de presión. Con el fin de abordar estas deficiencias, se introduce un nuevo conjunto de ecuaciones para mecánica del sólido formulada como un sistema de leyes de conservación de primer orden basada en una formulación mixta. Fundamentalmente, el nuevo sistema de leyes de conservación tiene una estructura similar a la de las famosas ecuaciones de Euler en el contexto de la Dinámica de Fluidos Computacional (CFD). Esto nos permite aprovechar algunas de las tecnologías CFD disponibles y adaptar el método en el contexto de la Mecánica de Sólidos. Esta tesis se basa en el trabajo realizado en Lee et al. 2013 mediante el desarrollo de la estructura de volúmenes finitos centrados en celdas upwind para el análisis numérico de dinámica del sólido explícita en grandes deformaciones y su implementación específicamente diseñada dentro del software de código abierto OpenFOAM, ampliamente utilizado ámbito académico e industrial. Además, la naturaleza orientada a objetos de su implementación proporciona una plataforma muy eficiente para su desarrollo posterior. En este marco computacional, las incógnitas básicas de este sistema son el momento lineal y el tensor gradiente de deformación. Asimismo, la formulación se extiende adicionalmente para un conjunto adicional de medidas de deformación que comprenden el cofactor del tensor gradiente de deformación y el jacobiano de deformación, con el fin de simular modelos constitutivos policonvexos que aseguran la estabilidad del material. El dominio se discretiza espacialmente usando un marco centrado en células de tipo Godunov estándar, donde se consigue la precisión de segundo orden empleando un procedimiento de reconstrucción lineal junto con un limitador de pendiente. Esto conduce a discontinuidades en las variables en la interfase de la célula que motivan el uso de un solucionador de Riemann mediante la introducción de un sesgo contra el viento en la evaluación de flujos de contacto numéricos. El presente solucionador acústico de Riemann es posteriormente desarrollado aplicando disipación pre-condicionada para mejorar su rendimiento en el cercano pero incompresibilidad régimen y extender su gama a aplicaciones de contacto. Además, se proponen dos marcos evolutivos en este estudio para satisfacer las involuciones subyacentes (o condiciones de compatibilidad) del sistema. Además, la discretización espacial se representa alternativamente a través de un marco de volumen finito centrado en células nodales para fines de comparación. Desde el punto de vista de la discretización temporal, se emplea un integrador temporal de Runge-Kutta de dos etapas con Disminución de Variación Total para asegurar segundo orden de precision. Finalmente, se simulan ejemplos numéricos de referencia para demostrar varios aspectos de la formulación que incluyen convergencia de malla, conservación de momento y la naturaleza libre de bloqueo de la formulación en dominios computacionales complejos.