Moving least‐squares particle hydrodynamics II: conservation and boundaries

Previous work by the author has shown that the consistency of the SPH method can be improved to acceptable levels by substituting MLS interpolants for SPH interpolants, that the SPH inconsistency drives the tension instability and that imposition of consistency via MLS severely retards tension instability growth. The new method however was not conservative, and made no provision for boundary conditions. Conservation is an essential property in simulations where large localized mass, momentum or energy transfer occurs such as high‐velocity impact or explosion modeling. A new locally conservative MLS variant of SPH that naturally incorporates realistic boundary conditions is described. In order to provide for the boundary fluxes one must identify the boundary particles. A new, purely geometric boundary detection technique for assemblies of spherical particles is described. A comparison with SPH on a ball‐and‐plate impact simulation shows qualitative improvement. Copyright © 2000 John Wiley & Sons, Ltd.     

Simulation of high velocity concrete fragmentation using SPH/MLSPH

The simulation of concrete fragmentation under explosive loading by a meshfree Lagrangian method, the smooth particle hydrodynamics method (SPH) is described. Two improvements regarding the completeness of the SPH‐method are examined, first a normalization developed by Johnson and Beissel (NSPH) and second a moving least square (MLS) approach as modified by Scheffer (MLSPH). The SPH‐Code is implemented in FORTRAN 90 and parallelized with MPI. A macroscopic constitutive law with isotropic damage for fracture and fragmentation for concrete is implemented in the SPH‐Code. It is shown that the SPH‐method is able to simulate the fracture and fragmentation of concrete slabs under contact detonation. The numerical results from the different SPH‐methods are compared with the data from tests. The good agreement between calculation and experiment suggests that the SPH‐program can predict the correct maximum pressure as well as the damage of the concrete slabs. Finally the fragment distributions of the tests and the numerical calculations are compared. Copyright © 2003 John Wiley & Sons, Ltd.     

Discontinuous Galerkin methods for non‐linear elasticity

This paper presents the formulation and a partial analysis of a class of discontinuous Galerkin methods for quasistatic non‐linear elasticity problems. These methods are endowed with several salient features. The equations that define the numerical scheme are the Euler–Lagrange equations of a one‐field variational principle, a trait that provides an elegant and simple derivation of the method. In consonance with general discontinuous Galerkin formulations, it is possible within this framework to choose different numerical fluxes. Numerical evidence suggests the absence of locking at near‐incompressible conditions in the finite deformations regime when piecewise linear elements are adopted. Finally, a conceivable surprising characteristic is that, as demonstrated with numerical examples, these methods provide a given accuracy level for a comparable, and often lower, computational cost than conforming formulations. Stabilization is occasionally needed for discontinuous Galerkin methods in linear elliptic problems. In this paper we propose a sufficient condition for the stability of each linearized non‐linear elastic problem that naturally includes material and geometric parameters; the latter needed to account for buckling. We then prove that when a similar condition is satisfied by the discrete problem, the method provides stable linearized deformed configurations upon the addition of a standard stabilization term. We conclude by discussing the complexity of the implementation, and propose a computationally efficient approach that avoids looping over both elements and element faces. Several numerical examples are then presented in two and three dimensions that illustrate the performance of a selected discontinuous Galerkin method within the class. Copyright © 2006 John Wiley & Sons, Ltd.     

Moving least‐square interpolants in the hybrid particle method

The hybrid particle method (HPM) is a particle‐based method for the solution of high‐speed dynamic structural problems. In the current formulation of the HPM, a moving least‐squares (MLS) interpolant is used to compute the derivatives of stress and velocity components. Compared with the use of the MLS interpolant at interior particles, the boundary particles require two additional treatments in order to compute the derivatives accurately. These are the rotation of the local co‐ordinate system and the imposition of boundary constraints, respectively. In this paper, it is first shown that the derivatives found by the MLS interpolant based on a complete polynomial are indifferent to the orientation of the co‐ordinate system. Secondly, it is shown that imposing boundary constraints is equivalent to employing ghost particles with proper values assigned at these particles. The latter can further be viewed as placing the boundary particle in the centre of a neighbourhood that is formed jointly by the original neighbouring particles and the ghost particles. The benefit of providing a symmetric or a full circle of neighbouring points is revealed by examining the error terms generated in approximating the derivatives of a Taylor polynomial by using a linear‐polynomial‐based MLS interpolant. Symmetric boundaries have mostly been treated by using ghost particles in various versions of the available particle methods that are based on the strong form of the conservation equations. In light of the equivalence of the respective treatments of imposing boundary constraints and adding ghost particles, an alternative treatment for symmetry boundaries is proposed that involves imposing only the symmetry boundary constraints for the HPM. Numerical results are presented to demonstrate the validity of the proposed approach for symmetric boundaries in an axisymmetric impact problem. Copyright © 2005 John Wiley & Sons, Ltd.     

A Lagrangian gradient smoothing method for solid‐flow problems using simplicial mesh

A novel Lagrangian gradient smoothing method (L‐GSM) is developed to solve “solid‐flow” (flow media with material strength) problems governed by Lagrangian form of Navier‐Stokes equations. It is a particle‐like method, similar to the smoothed particle hydrodynamics (SPH) method but without the so‐called tensile instability that exists in the SPH since its birth. The L‐GSM uses gradient smoothing technique to approximate the gradient of the field variables, based on the standard GSM that was found working well with Euler grids for general fluids. The Delaunay triangulation algorithm is adopted to update the connectivity of the particles, so that supporting neighboring particles can be determined for accurate gradient approximations. Special techniques are also devised for treatments of 3 types of boundaries: no‐slip solid boundary, free‐surface boundary, and periodical boundary. An advanced GSM operation for better consistency condition is then developed. Tensile stability condition of L‐GSM is investigated through the von Neumann stability analysis as well as numerical tests. The proposed L‐GSM is validated by using benchmarking examples of incompressible flows, including the Couette flow, Poiseuille flow, and 2D shear‐driven cavity. It is then applied to solve a practical problem of solid flows: the natural failure process of soil and the resultant soil flows. The numerical results are compared with theoretical solutions, experimental data, and other numerical results by SPH and FDM to evaluate further L‐GSM performance. It shows that the L‐GSM scheme can give a very accurate result for all these examples. Both the theoretical analysis and the numerical testing results demonstrate that the proposed L‐GSM approach restores first‐order accuracy unconditionally and does not suffer from the tensile instability. It is also shown that the L‐GSM is much more computational efficient compared with SPH, especially when a large number of particles are employed in simulation.     

A three‐dimensional element‐free Galerkin elastic and elastoplastic formulation

A small strain, three‐dimensional, elastic and elastoplastic Element‐Free Galerkin (EFG) formulation is developed. Singular weight functions are utilized in the Moving‐Least‐Squares (MLS) determination of shape functions and shape function derivatives allowing accurate, direct nodal imposition of essential boundary conditions. A variable domain of influence EFG method is introduced leading to increased efficiency in computing the MLS shape functions and their derivatives. The elastoplastic formulations are based on the consistent tangent operator approach and closely follow the incremental formulations for non‐linear analysis using finite elements. Several linear elastic and small strain elastoplastic numerical examples are presented to verify the accuracy of the numerical formulations. Copyright © 1999 John Wiley & Sons, Ltd.     

Topology optimization of structures using meshless density variable approximants

This paper proposes a new structural topology optimization method using a dual‐level point‐wise density approximant and the meshless Galerkin weak‐forms, totally based on a set of arbitrarily scattered field nodes to discretize the design domain. The moving least squares (MLS) method is used to construct shape functions with compactly supported weight functions, to achieve meshless approximations of system state equations. The MLS shape function with the zero‐order consistency will degenerate to the well‐known ‘Shepard function’, while the MLS shape function with the first‐order consistency refers to the widely studied ‘MLS shape function’. The Shepard function is then applied to construct a physically meaningful dual‐level density approximant, because of its non‐negative and range‐restricted properties. First, in terms of the original set of nodal density variables, this study develops a nonlocal nodal density approximant with enhanced smoothness by incorporating the Shepard function into the problem formulation. The density at any node can be evaluated according to the density variables located inside the influence domain of the current node. Second, in the numerical implementation, we present a point‐wise density interpolant via the Shepard function method. The density of any computational point is determined by the surrounding nodal densities within the influence domain of the concerned point. According to a set of generic design variables scattered at field nodes, an alternative solid isotropic material with penalization model is thus established through the proposed dual‐level density approximant. The Lagrangian multiplier method is included to enforce the essential boundary conditions because of the lack of the Kronecker delta function property of MLS meshless shape functions. Two benchmark numerical examples are employed to demonstrate the effectiveness of the proposed method, in particular its applicability in eliminating numerical instabilities. Copyright © 2012 John Wiley & Sons, Ltd.     

On elimination of shear locking in the element‐free Galerkin method

In this study, a method for completely eliminating the presence of transverse shear locking in the application of the element‐free Galerkin method (EFGM) to shear‐deformable beams and plates is presented. The matching approximation fields concept of Donning and Liu has shown that shear locking effects may be prevented if the approximate rotation fields are constructed with the innate ability to match the approximate slope (first derivative of displacement) fields and is adopted. Implementation of the matching fields concept requires the computation of the second derivative of the shape functions. Thus, the shape functions for displacement fields, and therefore the moving least‐squares (MLS) weight function, must be at least C¹ continuous. Additionally, the MLS weight functions must be chosen such that successive derivatives of the MLS shape function have the ability to exactly reproduce the functions from which they were derived. To satisfy these requirements, the quartic spline weight function possessing C² continuity is used in this study. To our knowledge, this work is the first attempt to address the root cause of shear locking phenomenon within the framework of the element‐free Galerkin method. Several numerical examples confirm that bending analyses of thick and thin beams and plates, based on the matching approximation fields concept, do not exhibit shear locking and provide a high degree of accuracy for both displacement and stress fields. Copyright © 2001 John Wiley & Sons, Ltd.     

Computational schemes for non‐linear elasto‐dynamics

This paper deals with the development of computational schemes for the dynamic analysis of non‐linear elastic systems. The focus of the investigation is on the derivation of unconditionally stable time‐integration schemes presenting high‐frequency numerical dissipation for these types of problem. At first, schemes based on Galerkin and time‐discontinuous Galerkin approximations applied to the equations of motion written in the symmetric hyperbolic form are proposed. Though useful, these schemes require casting the equations of motion in the symmetric hyperbolic form, which is not always possible. Furthermore, this approaches to unacceptably high computational costs. Next, unconditionally stable schemes are proposed that do not rely on the symmetric hyperbolic form. Both energy‐preserving and energy‐decaying schemes are derived. Numerical examples are presented to demonstrate the accuracy and efficiency of the proposed schemes. Copyright © 1999 John Wiley & Sons, Ltd.     

The SPH equations for fluids

Existing smoothed particle hydrodynamics (SPH) formulations for simulating continuous fluids have errors that may be divergent and it has been known for some time that the SPH equations do not satisfy low‐order polynomial completeness conditions. Here SPH equations are derived that have convergent error terms and a correction method is presented for enforcing low‐order polynomial completeness irrespective of how many completeness conditions are required. Discretization is achieved through division of the model domain, in its initial state, into sub‐domains that have Lagrangian boundaries. It is shown that boundary integrals appearing in one derivation of the SPH equations may be treated as a convergent error. In simulations of basic fluid flows convergence and zeroth‐order completeness are demonstrated, but significant instabilities and a failure to conserve energy are observed. Copyright © 2009 John Wiley & Sons, Ltd.     

A hybrid variational‐collocation immersed method for fluid‐structure interaction using unstructured T‐splines

We present a hybrid variational‐collocation, immersed, and fully‐implicit formulation for fluid‐structure interaction (FSI) using unstructured T‐splines. In our immersed methodology, we define an Eulerian mesh on the whole computational domain and a Lagrangian mesh on the solid domain, which moves arbitrarily on top of the Eulerian mesh. Mathematically, the problem reduces to solving three equations, namely, the linear momentum balance, mass conservation, and a condition of kinematic compatibility between the Lagrangian displacement and the Eulerian velocity. We use a weighted residual approach for the linear momentum and mass conservation equations, but we discretize directly the strong form of the kinematic relation, deriving a hybrid variational‐collocation method. We use T‐splines for both the spatial discretization and the information transfer between the Eulerian mesh and the Lagrangian mesh. T‐splines offer us two main advantages against non‐uniform rational B‐splines: they can be locally refined and they are unstructured. The generalized‐α method is used for the time discretization. We validate our formulation with a common FSI benchmark problem achieving excellent agreement with the theoretical solution. An example involving a partially immersed solid is also solved. The numerical examples show how the use of T‐junctions and extraordinary nodes results in an accurate, efficient, and flexible method. Copyright © 2015 John Wiley & Sons, Ltd.     

A novel non‐linearly explicit second‐order accurate L‐stable methodology for finite deformation: hypoelastic/hypoelasto‐plastic structural dynamics problems

A novel non‐linearly explicit second‐order accurate L‐stable computational methodology for integrating the non‐linear equations of motion without non‐linear iterations during each time step, and the underlying implementation procedure is described. Emphasis is placed on illustrative non‐linear structural dynamics problems employing both total/updated Lagrangian formulations to handle finite deformation hypoelasticity/hypoelasto‐plasticity models in conjunction with a new explicit exact integration procedure for a particular rate form constitutive equation. Illustrative numerical examples are shown to demonstrate the robustness of the overall developments for non‐linear structural dynamics applications. Copyright © 2004 John Wiley & Sons, Ltd.     

A refined non‐linear non‐conforming triangular plate/shell element

A refined non‐conforming triangular plate/shell element for geometric non‐linear analysis of plates/shells using the total Lagrangian/updated Lagrangian approach is constructed in this paper based on the refined non‐conforming element method for geometric non‐linear analysis. The Allman's triangular plane element with vertex degrees of freedom and the refined triangular plate‐bending element RT9 are used to construct the present element. Numerical examples demonstrate that the accuracy of the new element is quite high in the geometric non‐linear analysis of plates/shells. Copyright © 2003 John Wiley & Sons, Ltd.     

On the stability and convergence of a Galerkin reduced order model (ROM) of compressible flow with solid wall and far‐field boundary treatment

A reduced order model (ROM) based on the proper orthogonal decomposition (POD)/Galerkin projection method is proposed as an alternative discretization of the linearized compressible Euler equations. It is shown that the numerical stability of the ROM is intimately tied to the choice of inner product used to define the Galerkin projection. For the linearized compressible Euler equations, a symmetry transformation motivates the construction of a weighted L² inner product that guarantees certain stability bounds satisfied by the ROM. Sufficient conditions for well‐posedness and stability of the present Galerkin projection method applied to a general linear hyperbolic initial boundary value problem (IBVP) are stated and proven. Well‐posed and stable far‐field and solid wall boundary conditions are formulated for the linearized compressible Euler ROM using these more general results. A convergence analysis employing a stable penalty‐like formulation of the boundary conditions reveals that the ROM solution converges to the exact solution with refinement of both the numerical solution used to generate the ROM and of the POD basis. An a priori error estimate for the computed ROM solution is derived, and examined using a numerical test case. Published in 2010 by John Wiley & Sons, Ltd.     

An element‐wise,locally conservative Galerkin (LCG) method for solving diffusion and convection–diffusion problems

An element‐wise locally conservative Galerkin (LCG) method is employed to solve the conservation equations of diffusion and convection–diffusion. This approach allows the system of simultaneous equations to be solved over each element. Thus, the traditional assembly of elemental contributions into a global matrix system is avoided. This simplifies the calculation procedure over the standard global (continuous) Galerkin method, in addition to explicitly establishing element‐wise flux conservation. In the LCG method, elements are treated as sub‐domains with weakly imposed Neumann boundary conditions. The LCG method obtains a continuous and unique nodal solution from the surrounding element contributions via averaging. It is also shown in this paper that the proposed LCG method is identical to the standard global Galerkin (GG) method, at both steady and unsteady states, for an inside node. Thus, the method, has all the advantages of the standard GG method while explicitly conserving fluxes over each element. Several problems of diffusion and convection–diffusion are solved on both structured and unstructured grids to demonstrate the accuracy and robustness of the LCG method. Both linear and quadratic elements are used in the calculations. For convection‐dominated problems, Petrov–Galerkin weighting and high‐order characteristic‐based temporal schemes have been implemented into the LCG formulation. Copyright © 2007 John Wiley & Sons, Ltd.     

Non‐linear version of stabilized conforming nodal integration for Galerkin mesh‐free methods

A stabilized conforming (SC) nodal integration, which meets the integration constraint in the Galerkin mesh‐free approximation, is generalized for non‐linear problems. Using a Lagrangian discretization, the integration constraints for SC nodal integration are imposed in the undeformed configuration. This is accomplished by introducing a Lagrangian strain smoothing to the deformation gradient, and by performing a nodal integration in the undeformed configuration. The proposed method is independent to the path dependency of the materials. An assumed strain method is employed to formulate the discrete equilibrium equations, and the smoothed deformation gradient serves as the stabilization mechanism in the nodally integrated variational equation. Eigenvalue analysis demonstrated that the proposed strain smoothing provides a stabilization to the nodally integrated discrete equations. By employing Lagrangian shape functions, the computation of smoothed gradient matrix for deformation gradient is only necessary in the initial stage, and it can be stored and reused in the subsequent load steps. A significant gain in computational efficiency is achieved, as well as enhanced accuracy, in comparison with the mesh‐free solution using Gauss integration. The performance of the proposed method is shown to be quite robust in dealing with non‐uniform discretization. Copyright © 2002 John Wiley & Sons, Ltd.     

Weight‐adjusted discontinuous Galerkin methods: Matrix‐valued weights and elastic wave propagation in heterogeneous media

Weight‐adjusted inner products are easily invertible approximations to weighted L² inner products. These approximations can be paired with a discontinuous Galerkin (DG) discretization to produce a time‐domain method for wave propagation which is low storage, energy stable, and high‐order accurate for arbitrary heterogeneous media and curvilinear meshes. In this work, we extend weight‐adjusted DG methods to the case of matrix‐valued weights, with the linear elastic wave equation as an application. We present a DG formulation of the symmetric form of the linear elastic wave equation, with upwind‐like dissipation incorporated through simple penalty fluxes. A semidiscrete convergence analysis is given, and numerical results confirm the stability and high‐order accuracy of weight‐adjusted DG for several problems in elastic wave propagation.     

Smooth,second order,non‐negative meshfree approximants selected by maximum entropy

We present a family of approximation schemes, which we refer to as second‐order maximum‐entropy (max‐ent) approximation schemes, that extends the first‐order local max‐ent approximation schemes to second‐order consistency. This method retains the fundamental properties of first‐order max‐ent schemes, namely the shape functions are smooth, non‐negative, and satisfy a weak Kronecker‐delta property at the boundary. This last property makes the imposition of essential boundary conditions in the numerical solution of partial differential equations trivial. The evaluation of the shape functions is not explicit, but it is very efficient and robust. To our knowledge, the proposed method is the first higher‐order scheme for function approximation from unstructured data in arbitrary dimensions with non‐negative shape functions. As a consequence, the approximants exhibit variation diminishing properties, as well as an excellent behavior in structural vibrations problems as compared with the Lagrange finite elements, MLS‐based meshfree methods and even B‐Spline approximations, as shown through numerical experiments. When compared with usual MLS‐based second‐order meshfree methods, the shape functions presented here are much easier to integrate in a Galerkin approach, as illustrated by the standard benchmark problems. Copyright © 2009 John Wiley & Sons, Ltd.     

Higher‐order accurate time‐step‐integration algorithms by post‐integration techniques

In this paper, a framework to construct higher‐order‐accurate time‐step‐integration algorithms based on the post‐integration techniques is presented. The prescribed initial conditions are naturally incorporated in the formulations and can be strongly or weakly enforced. The algorithmic parameters are chosen such that unconditionally A‐stable higher‐order‐accurate time‐step‐integration algorithms with controllable numerical dissipation can be constructed for linear problems. Besides, it is shown that the order of accuracy for non‐linear problems is maintained through the relationship between the present formulation and the Runge–Kutta method. The second‐order differential equations are also considered. Numerical examples are given to illustrate the validity of the present formulation. Copyright © 2001 John Wiley & Sons, Ltd.