A Lagrangian particle method for the simulation of linear and nonlinear elastic models of soft tissue

doi:10.1016/j.jcp.2008.05.016

Journal of Computational Physics

Volume 227, Issue 21, 10 November 2008, Pages 9195-9215

https://doi.org/10.1016/j.jcp.2008.05.016 Get rights and content

Abstract

We present a novel Lagrangian particle method for the simulation of linear and nonlinear elastic models of soft tissue. Linear solids are represented by the Lagrangian formulation of the stress–strain relationship that is extended to nonlinear solids by using the Lagrangian evolution of the deformation gradient described in a moving framework. The present method introduces a level set description, along with the particles, to capture the body deformations and to enforce the boundary conditions. Furthermore, the accuracy of the method in cases of large deformations is ensured by implementing a particle remeshing procedure. The method is validated in several benchmark problems, in two and three dimensions and the results compare well with the results of respective finite elements simulations. In simulations of large solid deformation under plane strain compression, the finite element solver exhibits spurious structures that are not present in the Lagrangian particle simulations. The particle simulations are compared with experimental results in an aspiration test of liver tissue.

Introduction

The simulation of soft tissue is one of the key components of virtual reality systems with applications ranging from video games to virtual surgery environments [48], [47]. An important aspect of these simulations is the physics-based modeling of linear and nonlinear elastic solids undergoing large deformations. A number of computational techniques have been employed in the past in order to address this problem including finite element [7], [40], [31], finite differences [49], mass–spring models [50], [51] and particle methods [15], [16], [33], [38], [39]. Grid based methods such as finite elements have been shown to be efficient and robust in simulations of systems undergoing small or medium deformations [7], whereas meshless/particle methods are advocated for the simulations of solids undergoing excessive deformation and mechanical splitting [38].

Meshless methods represent material properties on particles and the evolution of the material is linked to the evolution of the computational elements. They are usually based on a Galerkin or collocation discretization of the governing equations and may rely only on particles (as in the case of the Reproducing Kernel Particle Method (RKPM) [30]) or they may in addition employ a grid for the computation of the stress–strain relations (Material Point Method (MPM) [45]). The adaptivity and grid-free character of meshless methods have enabled a number of state of the art simulation in biomechanics [9]. We note however that, even though meshless methods are considered capable of handling large deformations, they often rely on empirical techniques in order to maintain the accuracy of the method. This problem may be further aggravated when convection terms distort the computational elements, as in simulations of fluid flows, that may be viewed as a mechanical system undergoing extreme deformations. Smoothed Particle Hydrodynamics (SPH) has been proposed as a meshless method capable of resolving problems with large deformations for both solids and fluids. In SPH the continuum properties are discretized on smooth particles, the stress–strain governing equations are formulated in a Lagrangian frame and the derivatives are computed by taking the derivatives of the particle kernels. SPH allows for a unifying formulation of stress–strain relations and has been used extensively in order to simulate complex systems ranging from astrophysics [13], [32] to impact simulations in solid mechanics [38]. The simulations of solids using SPH suffer from the well-known tensile instability. The problem of tensile instability arises when the distance of particles is small under positive pressure, an effect attributed to the shape of the second derivative of the interpolating kernel [46]. The forces become attractive due to the shape of the derivative approximation resulting in large numerical errors. A number of techniques have been presented to resolve this problem including the use of an alternative smoothing kernel [24], the introduction of stress points [39] or artificial stress [33], [15]. We consider that this pathology of SPH is related to the fact that smooth particle approximations are inaccurate when the particle distribution is excessively distorted [2], [5]. In fluids problems convection is often responsible for large particle distortions whereas in solid mechanics particle distortion is associated with a nonlinear strain field. Several works [38], [15], [33] have employed particle models for linear elasticity but the use of particle methods in problems with nonlinear elasticity is relatively limited.

Nonlinear elasticity problems have been simulated in an Eulerian framework using the Finite Element Method (FEM) [40]. These simulations are based on the formulation of the deformation gradient in terms of the current material position with respect to its reference position. Similarly, the Reproducing Kernel Particle Method (RKPM) solve the governing equations of nonlinear elasticity in the framework of the undeformed (reference) configuration [6]. The extension of this approach to Lagrangian particle methods is hindered by the fact that when the material properties are tracked in a moving framework the reference position becomes irrelevant. In order to circumvent this difficulty we replace the direct evaluation of the deformation gradient by its temporal evolution. The resulting formulation enables a straightforward and efficient implementation using particle methods. To the best of our knowledge, particle models for nonlinear elasticity have not been presented in the literature. In fluid simulations, the nonlinear relationship between stress tensor and velocity gradient has been studied recently by the use of non-Newtonian models such as Bingham, Bulkley, and power law models [22], [28].

The proposed particle method, relies on remeshing to ensure the accuracy of the simulations in particular for solids undergoing large deformations. Particles are convected in a Lagrangian framework, followed by a regularization of their locations and the corresponding projection of the particle properties. This approach enables accurate simulations of solids undergoing large deformations and eliminates spurious effects such as the tensile instability. In addition, this remeshing procedure accommodates the generalisation of stress–strain relations in a Lagrangian framework and the extension of particle methods to problems of nonlinear elasticity. We also introduce a particle-level set description [19], along with the particles, to capture the solid bodies and to enforce the boundary conditions.

We compare the present particle methods with the finite element solver (Abaqus 6.4/Explicit) in benchmark problems involving linear and nonlinear elastic solids. We examine the accuracy of the method in a three-dimensional test problem with periodic boundary conditions and demonstrate its capability to simulate large deformations in a plane strain compression test. Finally, we validate our particle solver by simulating liver tissue exposed to an aspiration test [35], [36].

The paper is structured as follows: in Section 2 we introduce the governing equations for linear and nonlinear elastic material followed by their particle representation in Sections 3 Particle representation of solids, 4 Particle equations. We validate the particle solver in comparison to an analytical solution (Section 5) and in a benchmark problem along with the FEM solution (Section 6) for both, linear and nonlinear elasticity. In Section 7 we present the simulation of the mechanical behavior of liver tissue during an aspiration test and compare with experimental data. Moreover, we consider the simulation problem of a rotating elastic cylinder (Section 8) and of a bouncing ball (Section 9). The paper concludes with an assessment of the method and outlines future work.

Section snippets

Governing equations

We consider the Lagrangian formulation of the governing equations for a linear and a nonlinear model of an elastic material.

Function and gradient approximations using particles

In the context of particle methods [17], [5], [27] a smooth approximation of a function $Φ (x)$ can be constructed by using a mollification kernel $ζ_{ϵ} (x)$ : $Φ_{ϵ} (x) = Φ ★ ζ_{ϵ} = \int Φ (y) ζ_{ϵ} (x - y) d y,$ where $ϵ$ denotes a characteristic length of the kernel.

The kernel is of order r when the following moment conditions are satisfied: $\int ζ_{ϵ} (x) d x = 1,$ $\int x^{i} ζ_{ϵ} (x) d x = 0 if | i | ⩽ r - 1,$ $\int | x |^{r} ζ_{ϵ} (x) d x ⩽ \infty .$ This mollified approximation $Φ_{ϵ} (x)$ can be discretized using the particle locations as quadrature points and a particle approximation of the

Particle equations

The particle position $x_{p}$ , mass $m_{p}$ , volume $v_{p}$ , and velocity $u_{p}$ evolve by the following system of ordinary differential equations: $\begin{matrix} \frac{d x_{p}}{d t} = u_{p}, \\ \frac{d m_{p}}{d t} = 0, \\ \frac{d v_{p}}{d t} = 〈 \nabla \cdot u 〉_{p} v_{p}, \end{matrix}$ $\frac{d u_{p}}{d t} = \frac{v_{p}}{m_{p}} (- 〈 \nabla p 〉_{p} + 〈 \nabla \cdot S 〉_{p}),$ where $〈 \nabla ♢ 〉_{p}$ represents the mollified derivative approximation of $♢$ based on Eqs. (29), (30). The evaluation of the pressure p and the deviatoric stress S depends on the constitutive model of the elastic material.

The surface of the elastic body is described using the Particle Level Set Method [18], [19].

The Level

Accuracy

We consider an elastic unit cube with periodic boundary conditions in all directions. A sinusoidal velocity $(u (x, 0) = \sin (2 π x) \sin (2 π y) \sin (2 π z))$ is initially imposed that leads to periodic oscillations in the material. We quantify the accuracy of the particle method by the maximum error $(L_{\infty})$ of the velocity profile as compared with the analytic solution after one period. Young’s modulus is $E = 100$ and the Poisson ratio is set to $ν = 0.3$ in the linear case. The polynomial coefficient is $C_{10} = 10$ and the

Plane strain compression test

We compare the results of plane strain compression test with an FEM solution provided by Abaqus 6.4/Explicit for linear and nonlinear elasticity. The elastic material is initially undeformed and has a rectangular shape of size $1 \times 2$ (Fig. 2). The horizontal faces move with a constant vertical velocity $u_{piston} = 25$ whereas the vertical faces are considered as free surfaces. The vertical compression of the material leads to a horizontal expansion. The thickness of the piston is neglected. The

Simulation of aspiration test on liver tissue

The particle method is applied to the simulation of the deformation of liver tissue during an aspiration test [35], [36]. The aspiration test is performed by pressing a tube against the liver tissue, creating a vacuum so that the tissue is sucked into the aspiration cavity (Fig. 6, Fig. 7). The material properties of the tissue are determined by measuring the tissue height inside the tube.

We adopt the hyperelastic liver model of Nava et al. [35] that describes the strain–energy relationship in

Simulation of a rotating elastic cylinder

The simulation of a rotating elastic cylinder allows for an assessment of the conservation of angular momentum in SPH, an issue that has been thoroughly addressed in a work of Hoover et al. [21]. In the present work the linear elastic cylinder of radius $r_{0} = 1$ is initially undeformed and rotates around its symmetry axis with the angular speed $ω = 2 π$ . The material of the cylinder has Young’s modulus of $E = 100$ and Poisson ratio of $ν = 0.3$ . The time integration involves the Runge–Kutta scheme of 4th

Simulation of a bouncing ball

A linear elastic spherical ball of radius $r = 0.28$ hits a fixed wall at $t = 0$ . The ball has an initial velocity of $u_{0} (x) = 1$ and Young’s modulus $E = 100$ and Poisson ratio $ν = 0.4$ . The time integration scheme is the Runge–Kutta scheme of 4th order with a time step $Δ t = 1 \times 10^{- 4}$ .

The particles are remeshed every 10th time step, and the level sets are reinitialized every 500 time steps. Fig. 14 shows the deformation of the ball discretized by 8000 particles at $t = 0.00, 0.03, 0.07, 0.14, 0.21$ and 0.25, while Fig. 15

Conclusions

In this paper we present a novel particle method for the simulation of linear and nonlinear elastic solid models. The particle method relies on a remeshing procedure to ensure convergence and to avoid numerical artifacts such as tensile instabilities that are often associated with SPH formulations.

The formulation enables simulations of linear and nonlinear elastic materials. In the latter case we introduce, to the best of our knowledge, the first ever Lagrangian particle simulations of

Acknowledgments

We would like to thank the group of Prof. Mazza (IMES, ETHZ) for several helpful discussions. This project was funded by the Swiss National Science Foundation, NCCR ‘Computer Aided and Image Guided Medical Interventions’ (Co-Me).

The paper is dedicated to the Tony Leonard for his guidance over the years and his contributions to our work. He has been a source of inspiration in every aspect of our research activities.

References (51)

A.K. Chaniotis et al.
Remeshed smoothed particle hydrodynamics for the simulation of viscous and heat conducting flows
J. Comput. Phys.
(2002)
G.-H. Cottet
Artificial viscosity models for vortex and particle methods
J. Comput. Phys.
(1996)
M. Doblare et al.
On the employ of meshless methods in biomechanics
Comput. Meth. Appl. Mech. Eng.
(2005)
Jeff D. Eldredge et al.
A vortex particle method for two-dimensional compressible flow
J. Comput. Phys.
(2002)
Jeff D. Eldredge et al.
A general deterministic treatment of derivatives in particle methods
J. Comput. Phys.
(2002)
J.P. Gray et al.
SPH elastic dynamics
Comput. Meth. Appl. Mech. Eng.
(2001)
Simone Elke Hieber et al.
A Lagrangian particle level set method
J. Comput. Phys.
(2005)
A. Khayyer et al.
Corrected incompressible SPH method for accurate water-surface tracking in breaking waves
Coastel Eng.
(2008)
J.J. Monaghan
SPH without a tensile instability
J. Comput. Phys.
(2000)
Stanley Osher et al.
Front propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulation
J. Comput. Phys.
(1988)

P.W. Randles et al.

Smoothed particle hydrodynamics: some recent improvements and applications

Comput. Meth. Appl. Mech. Engng.

(1996)

I.F. Sbalzarini et al.

PPM – a highly efficient parallel particle-mesh library for the simulation of continuum systems

J. Comput. Phys.

(2006)

Deborah Sulsky et al.

Application of a particle-in-cell method to solid mechanics

Comput. Phys. Commun.

(1995)

J.W. Swegle et al.

Smoothed particle hydrodynamics stability analysis

J. Comput. Phys.

(1995)

J.T. Beale

A convergent 3-D vortex method with grid-free stretching

Math. Comput.

(1986)

R.M. Christensen

Theory of Viscoelasticity

(2003)

G.-H. Cottet et al.

Vortex Methods – Theory and Practice

(2000)

K.T. Danielson et al.

Parallel computation of meshless methods for explicit dynamic analysis

Int. J. Numer. Methods Eng.

(2000)

G. Debunne, M. Desbrun, M.-P. Cani, A. Barr, Dynamic real-time deformations using space and timing adaptive sampling,...

P. Degond et al.

The weighted particle method for convection-diffusion equations. Part 1: The case of an isotropic viscosity

Math. Comput.

(1989)

Y.C. Fung

Biomechanics, Mechanical Properties of Living Tissues

(1993)

R.A. Gingold et al.

Smoothed particle hydrodynamics: theory and application to non-spherical stars

Mon. Not. R. Astron. Soc.

(1977)

R.A. Gingold et al.

The reliability of finite-difference and particle methods for fragmentation problems

Mon. Not. R. Astron. Soc.

(1982)

J.E. Guilkey et al.

Computational modeling of multicellular constructs with the material point method

J. Biomech.

(2006)

Ole H. Hald

Convergence of vortex methods for Euler’s equations, III

SIAM J. Numer. Anal.

(1987)

Cited by (28)

A remeshed vortex method for mixed rigid/soft body fluid–structure interaction
2021, Journal of Computational Physics
We outline a 2D algorithm for solving incompressible flow–structure interaction problems for mixed rigid/soft body representations, within a consistent framework based on the remeshed vortex method. We adopt the one-continuum formulation to represent both solid and fluid phases on an Eulerian grid, separated by a diffuse interface. Rigid solids are treated using Brinkman penalization while an inverse map technique is used to obtain elastic stresses in the hyperelastic solid phase. We test our solver against a number of benchmark problems, which demonstrate physical accuracy and first to second order convergence in space and time. Benchmarks are complemented by additional investigations that illustrate the ability of our numerical scheme to capture essential fluid–structure interaction phenomena across a variety of scenarios involving internal muscular actuation, self propulsion, multi-body contact, heat transfer and rectified viscous streaming effects. Through these illustrations, we showcase the ability of our solver to robustly deal with different constitutive models and boundary conditions, solve disparate multi-physics problems and achieve faster time-to-solutions by sidestepping CFL time step restrictions.
Thermal simulation in multiphase incompressible flows using coupled meshfree and particle level set methods
2018, Computer Methods in Applied Mechanics and Engineering
Citation Excerpt :
In accordance with the first utilization of remeshing in PLS methods [41], a brief description of this procedure is herein revisited. Although remeshing can directly be applied to physical particles and have successfully been used in some related works (e.g., [67,70]), this approach may however cause some difficulties in handling wall boundary conditions which would require additional considerations. This matter is thus not taken into account for the present work.
A particle-based numerical solver is presented, applicable to the simulation of heat transfer in multiphase immiscible flows including surface tension. In the context of meshfree methods, the Laplacian operator is recognized as the most numerically challenging ingredient of the heat equation. The well-known difficulty of approximating higher-order spatial derivatives with meshfree methods is herein addressed by adopting two advanced schemes in order to ensure second-order accuracy. In addition, a Lagrangian particle level set method with second-order reinitializationis introduced, for the first time in thermal simulations, to capture the location of the interface, i.e., moving and/or deformable boundaries of the continuum at each time-step. This leads to a more conservative solution, alleviating the mass loss during the simulation. Furthermore, a narrow band of the exterior geometric particles is exploited to form physical ghost particles to enforce the required boundary conditions. This novel approach ensures solver performance without the necessity of defining extra dummy particles for treating boundary conditions in meshfree simulations. Three benchmark problems are considered for evaluating the performance of the proposed solver against a reference analysis performed in COMSOL Multiphysics ^®.
Modeling extracellular matrix viscoelasticity using smoothed particle hydrodynamics with improved boundary treatment
2017, Computer Methods in Applied Mechanics and Engineering
Citation Excerpt :
A similar implementation using ALE or the immersed boundary method has not been published yet on migration of a cell through an ECM which requires solid–solid interaction and degradation of the ECM. Meshless particle based methods can be a valuable alternative here because they can deal with deformable interfaces [10], large deformations and discontinuities, while at the same time allowing a natural and efficient coupling with other particle methods such as the Discrete Element Method which can realistically model deformable cells as individual, discrete objects [11–13]. To the best of our knowledge, no computational model has been published that captures cell migration through a continuous ECM that functions as a mechanical obstacle for the cell and requires degradation for the cell to pass through.
We propose viscoelastic smoothed particle hydrodynamics (SPH) with extended boundary conditions as a new method to model the extracellular matrix (ECM) in contact with a migrating cell. The contact mechanics between a cell and ECM is modeled based on an existing boundary method in SPH that corrects for the well-known missing kernel support problem in Fluid Structure Interactions (FSI). This boundary method is here extended to allow the modeling of moving boundaries in contact with a viscoelastic solid. To validate the method, simulations are performed of tractions applied to a viscoelastic solid, Stokes flow around an array of square pillars, and indentation of a viscoelastic material with a circular indenter. By drop out of the inertial terms in the SPH equations of motion, the new SPH formulation allows to solve problems in a low Reynolds environment with a timestep independent of the particle spacing, permitting to model processes at the cellular scale (i.e. $μ m$ -scale). The potential of the method to capture cell–ECM interactions is demonstrated by simulation of a self propelling object that locally degrades the ECM by fluidizing it to permit migration. This should enable us to model and understand realistic cell–matrix interactions in the future.
A comparative study of penalization and phase field methods for the solution of the diffusion equation in complex geometries
2015, Journal of Computational Physics
Citation Excerpt :
Particle methods [43–45] have been proposed for the solution of PDEs such as diffusion equations, starting with the method of Particle Strength Exchange [46]. The incorporation of BC in such methods remains a topic of active research [11,47–49] with significant success in the formalism of the Reproducing Kernel Particle methods [50]. In recent years, it has become evident that particle methods simulating mechanical systems with large deformations need to be coupled with remeshing procedures [51,52] that employ a Cartesian grid.
We present a comparative study of penalization and phase field methods for the solution of the diffusion equation in complex geometries embedded using simple Cartesian meshes. The two methods have been widely employed to solve partial differential equations in complex and moving geometries for applications ranging from solid and fluid mechanics to biology and geophysics. Their popularity is largely due to their discretization on Cartesian meshes thus avoiding the need to create body-fitted grids. At the same time, there are questions regarding their accuracy and it appears that the use of each one is confined by disciplinary boundaries. Here, we compare penalization and phase field methods to handle problems with Neumann and Robin boundary conditions. We discuss extensions for Dirichlet boundary conditions and in turn compare with methods that have been explicitly designed to handle Dirichlet boundary conditions. The accuracy of all methods is analyzed using one and two dimensional benchmark problems such as the flow induced by an oscillating wall and by a cylinder performing rotary oscillations. This comparative study provides information to decide which methods to consider for a given application and their incorporation in broader computational frameworks. We demonstrate that phase field methods are more accurate than penalization methods on problems with Neumann boundary conditions and we present an error analysis explaining this result.
Multiscale modeling in food engineering
2013, Journal of Food Engineering
Since many years food engineers have attempted to describe physical phenomena such as heat and mass transfer that occur in food during unit operations by means of mathematical models. Foods are hierarchically structured and have features that extend from the molecular scale to the food plant scale. In order to reduce computational complexity, food features at the fine scale are usually not modeled explicitly but incorporated through averaging procedures into models that operate at the coarse scale. As a consequence, detailed insight into the processes at the microscale is lost, and the coarse scale model parameters are apparent rather than physical parameters. As it is impractical to measure these parameters for the large number of foods that exist, the use of advanced mathematical models in the food industry is still limited. A new modeling paradigm – multiscale modeling – has appeared that may alleviate these problems. Multiscale models are essentially a hierarchy of sub-models which describe the material behavior at different spatial scales in such a way that the sub-models are interconnected. In this article we will introduce the underlying physical and computational concepts. We will give an overview of applications of multiscale modeling in food engineering, and discuss future prospects.
A generalized wall boundary condition for smoothed particle hydrodynamics
2012, Journal of Computational Physics
Citation Excerpt :
For this reason SPH has a high potential especially for simulating multi-phase systems and can be applied to a broad variety of problems. Over the past three decades, SPH was successfully used to simulate complex problems ranging from magnetohydrodynamics [3] and solid mechanics [4–7] to fluid mechanics including free surfaces [8,9], surface tension [10,11] and transport phenomena [12,13]. Regardless of the application, boundary conditions are one of the key aspects of a numerical simulation and special attention should be paid to a correct and accurate representation of them.
In this paper we present a new formulation of the boundary condition at static and moving solid walls in SPH simulations. Our general approach is both applicable to two and three dimensions and is very simple compared to previous wall boundary formulations. Based on a local force balance between wall and fluid particles we apply a pressure boundary condition on the solid particles to prevent wall penetration. This method can handle sharp corners and complex geometries as is demonstrated with several examples. A validation shows that we recover hydrostatic equilibrium conditions in a static tank, and a comparison of the classical dam break simulation with state-of-the-art results in literature shows good agreement. We simulate various problems such as the flow around a cylinder and the backward facing step at Re = 100 to demonstrate the general applicability of this new method.

View all citing articles on Scopus

View full text

A Lagrangian particle method for the simulation of linear and nonlinear elastic models of soft tissue

Abstract

Introduction

Section snippets

Governing equations

Function and gradient approximations using particles

Particle equations

Accuracy

Plane strain compression test

Simulation of aspiration test on liver tissue

Simulation of a rotating elastic cylinder

Simulation of a bouncing ball

Conclusions

Acknowledgments

J. Comput. Phys.

J. Comput. Phys.

Comput. Meth. Appl. Mech. Eng.

J. Comput. Phys.

J. Comput. Phys.

Comput. Meth. Appl. Mech. Eng.

J. Comput. Phys.

Coastel Eng.

J. Comput. Phys.

J. Comput. Phys.

Comput. Meth. Appl. Mech. Engng.

J. Comput. Phys.

Comput. Phys. Commun.

J. Comput. Phys.

A convergent 3-D vortex method with grid-free stretching

Math. Comput.

Theory of Viscoelasticity

Vortex Methods – Theory and Practice

Parallel computation of meshless methods for explicit dynamic analysis

Int. J. Numer. Methods Eng.

The weighted particle method for convection-diffusion equations. Part 1: The case of an isotropic viscosity

Math. Comput.

Biomechanics, Mechanical Properties of Living Tissues

Smoothed particle hydrodynamics: theory and application to non-spherical stars

Mon. Not. R. Astron. Soc.

The reliability of finite-difference and particle methods for fragmentation problems

Mon. Not. R. Astron. Soc.

Computational modeling of multicellular constructs with the material point method

J. Biomech.

Convergence of vortex methods for Euler’s equations, III

SIAM J. Numer. Anal.