A Fully Explicit Lattice Spring Elastoplastic Solid Model for Fluid Structure Interaction Applications

The Fluid Structure Interaction (FSI) solver based on mesh-less methods developed at University of Nottingham Malaysia has been proven to be successful in solving hydroelastic problems involving elastic body. In this paper, the possibility of modelling elastocplastic material behaviour using mesh-less method is further explored. Elastoplastic material bahaviour is essential for modelling ductile fracture in a wide range of engineering problems. The use of a type of nonlocal lattice spring model (LSM) based on fully explicit integration scheme in modelling elastoplastic material behaviour is presented in this work. Specifically, the return mapping algorithm is incorporated in the solid solver for modelling J2 plasticity. The effect of time step size and particle spacing on the accuracy of the solver are investigated by checking the predicted strain against theoretical solution. It is found out that the strain results obtained by varying the time step size are somewhat consistent with the theoretical solution. The predicted strain also converges to the theoretical result as the particle spacing is refined.


Introduction
Ductile fracture in metals is critical in engineering application.In contrast to brittle fracture, ductile fracture is accompanied with plastic deformation and reduction in material load carrying capacity.Thus, to numerically model ductile fracture, it is important for the numerical model to incorporate plastic deformation.Conventional continuum-based method such as finite element method (FEM) has been widely applied in engineering practice.However, the challenge for FEM in simulating fracture lies in the continuous requirement of the model.Additionally, remeshing, prior crack path knowledge are among the necessary tools in realizing the crack propagation in FEM [1].Another class of method known as discrete-based method such as LSM is more suitable in modelling fracture as fracture can be realized by bond breakage governed by preset failure criteria.LSM has been applied extensively in the study of brittle fracture, for e.g.[2,3].Furthermore, the use of lattice type model to model ductile fracture has been developed and studied in [4].Recently, a type of LSM known as volume compensated particle method (VCPM) utilizing nonlocal interaction apart from local axial interaction is extended to model plasticity and ductile fracture [5].
In the work [5], VCPM is solved using an implicit numerical scheme known as atomic finite element method.Here, the symplectic Euler scheme is used to solve the equation of motion and the effect of varying particle and time step size on the accuracy of VCPM in solving elastoplastic test case is investigated.Return mapping algorithm is implemented in VCPM for modelling the elastoplastic material behaviour.Details of the implemented return mapping algorithm in VCPM and the 3D tensile beam test case used for validation study are described next.

Methodology
In this work, 3D simple cubic (SC) lattice structure consisting of six first-nearest neighbours and twelve second-nearest neighbours is adopted in the simulation.Under homogeneous isotropic material assumption, the interaction force between two VCPM solid particles is [6,7]: ))  ̂ where   is the total potential energy of the unit cell,   is half of the length change in spring connecting particle i and j,  ̂ is the unit vector between particle i and j, k is the local spring stiffness and T is the nonlocal spring potential.
The implemented J2 plasticity in VCPM is based on stress tensor which can be calculated from the interaction force shown in equation ( 1) [8].In other words, the stress tensor at a solid particle is formulated based on the spring bond deformation.Return mapping or elastic predictor plastic corrector algorithm is implemented along with J2 plasticity and isotropic hardening in this study.J2 or maximum distortion energy criterion disregards the effect of hydrostatic stress components on the yielding of material.J2 plasticity theory is described followed by the implementation of the model in VCPM.
The deviatoric stress tensor in index notation can be represented as: where   is the deviatoric stress tensor,   is the Kronecker delta and  = 1 3 is the mean hydrostatic stress.Subsequently, the yield function is defined as: where is the second invariant of deviatoric stress tensor,  is the initial yield strength of the material,  is the hardening modulus and ̅  is the accumulated plastic strain.Associative plastic flow rule is used where the yield function is taken as the plastic potential: where  ̇ is the plastic multiplier and ̅ is the equivalent von Mises stress.The Kuhn-Tucker consistency condition describing the loading and unloading conditions is satisfied by the model: From the consistency condition, when the yield function is less than zero, i.e., elastic deformation, the plastic multiplier is zero whereas if there is plastic flow, the yield function is zero from the constraint.Following the additive decomposition of total strain in plasticity theory, the spring bond deformation in VCPM also consists of elastic and plastic/permanent deformation: At the elastic predictor stage, the trial bond force is calculated using the elastic bond deformation.The trial yield function is then calculated as:  (7) If the trial yield function exceeds the yield surface, the plastic multiplier is larger than zero and can be determined as [5]: where  is the shear modulus.With the plastic multiplier, the incremental plastic strain tensor is updated as: and also, the equivalent plastic strain is updated as: Finally, the incremental plastic stretch can be calculated from the incremental plastic strain tensor as: where   is the original bond length and note that summation convention is implied in equation ( 11) for capital letter I and J.The corrected elastic bond elongation can then be obtained and used to recalculate the bond force after the plastic deformation is updated.The time step size is governed by the speed of P-wave: where  is the Courant-Friedrich-Levy condition and  is the particle spacing.  = √ +(4/3)  is the sound speed in solid with  being the bulk modulus and  the density of the solid.

Results and Discussions
The geometric details of the elastoplastic test case considered is shown in figure 1.The left and back face of the beam are fixed in x-and y-direction respectively.The top face of the beam is fixed in zdirection and the bottom face is being applied uniformly distributed load.The material constants considered here are Young's modulus E=146 MPa, Poisson ratio v=0.3, hardening modulus H=60.6 MPa and yield strength of 200000 Pa.The loading scenario for the first case is shown in figure 2. First, the load increases linearly from 0 to 250000 Pa at t=2.5 s.Then, unloading occurs from 250000 Pa to 200000 Pa before reloading occurs until 300000 Pa.To check the effect of time step size on the accuracy of the elastoplastic solver, the CFL value in equation ( 12) is varied and the corresponding stress strain curve is shown in figure 3. The particle spacing Dp used is 1.0 mm.From figure 3, it can be seen that as the time step size decreases, the strain predicted from numerical solution does not differ much and is consistent with the theoretical solution.This shows that the time step size criterion is suitable and applicable for elastoplastic scenario.Next, the particle spacing is refined and the convergence of the elastoplastic solver is checked.For the second loading scenario, the load is monotonically increased up to 250000 MPa.The particle spacings considered in this second case are Dp = 1.0 mm and Dp = 0.5 mm.From the stress strain curve shown in figure 4, when the particle spacing is refined, the predicted curve converges to theoretical result.

Conclusions
Modelling of elastoplastic material bahaviour using a type of nonlocal LSM known as VCPM has been demonstrated.Return mapping algorithm along with explicit integration scheme are adopted for the VCPM elastoplastic solver.It has been shown through the 3D tensile beam test case that the effect of varying the time step size is minimal on the accuracy of the solver.Moreover, as the particle spacing is refined, the predicted numerical strain converges to theoretical solution.With the validated elastoplastic solver, future work may focus on ductile fracture modelling and its application on real world scenario.

Figure 1 .
Figure 1.Geometric details the beam subjected to uniformly distributed load.

Figure 2 .
Figure 2. Loading history of uniformly distributed load.

Figure 3 .
Figure 3. Stress strain curve for different CFL values and the theoretical result.

Figure 4 .
Figure 4. Comparison of stress strain curve for theoretical and VCPM result obtained using Dp=1.0 mm and Dp=0.5 mm.