Material point method after 25 years: Theory, implementation, and applications

Alban de Vaucorbeil , ... Jian Ying Wu , in Advances in Applied Mechanics, 2020

2.1.1 Motion and deformation

In continuum mechanics, a body B is considered as being formed by an infinite set of material points, which are endowed with certain mechanical properties. The position vector of a material point in the initial, undeformed configuration of the body is denoted X relative to some coordinate basis. X is named the material or Lagrangian coordinate. The position of the same material point, in the deformed configuration, is designated by x called spatial or Eulerian coordinate.

The motion (deformation) of a solid is described by a function ϕ (X, t). A relation between spatial coordinates and material coordinates can be established as follows

(4) x = ϕ ( X , t )

The displacement, velocity and acceleration fields of a body are the primary kinematical fields in describing the motion of the body. The displacement of a material point X, denoted by u(X, t), is the difference between its current position ϕ(X, t) and its initial position ϕ(X, 0). So,

(5) u ( X , t ) : = ϕ ( X , t ) ϕ ( X , 0 ) = x X

The velocity of a material point X, denoted by v(X, t), is defined as the rate of change of position of this material point, that is

(6) v ( X , t ) : = ϕ ( X , t ) t

This is the Lagrangian velocity field. There exists a Eulerian form for this velocity field but as the MPM adopts a Lagrangian description, it is not discussed here.

The acceleration of a material point X is the rate of change of its velocity, i.e., the material time derivative of the velocity,

(7) a ( X , t ) : = v ( X , t ) t = 2 ϕ ( X , t ) t 2

The deformation gradient tensor F is a key quantity in finite deformation continuum mechanics as all deformation quantities are derived from it. It is a linear mapping operator which maps each infinitesimal linear element dX in the reference configuration into an infinitesimal linear element dx in the current configuration. It is defined as:

(8) F : = ϕ X = x X or F i j = x i X j

Next, the concept of material time derivative is introduced. To understand this important concept, consider the following situation. Assume we have a certain field φ (scalar, vectorial or tensorial) defined over the body for which we want to know the rate of change, at a given material point X. This is known as the material time derivative of φ. There are two definition of this concept, corresponding to material and spatial descriptions, respectively:

1.

Lagrangian description. In the Lagrangian description, the independent variables are the material coordinates X and time t. So all we have to do is taking the partial derivative of the given field φ with respect to time. For a material field φ(X, t), its material time derivative is

(9) D φ ( X , t ) D t φ . = φ ( X , t ) t

where the first two equations indicate standard notation for the material time derivative. As the MPM adopts a Lagrangian description, the material time derivative is very simple.
2.

Eulerian description. The considered field is φ(x, t). This case is much more complicated since not only time changes but also the spatial position x of the considered particle. We must calculate the partial derivative with respect to time of the material description of φ(x, t), keeping X fixed.

(10) D φ ( x , t ) D t φ . : = lim Δ t 0 φ ( ϕ ( X , t + Δ t ) , t + Δ t ) φ ( ϕ ( X , t ) , t ) Δ t

Using the chain rule, we obtain the important formula of the material time derivative for a Eulerian scalar field:

(11) D φ ( x , t ) D t = φ ( x , t ) t + φ ( x , t ) v ( x , t )

The term ∂φ/∂t is called the spatial time derivative, and the term φ , j v j is the convective term, which is also called the transport term.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/S0065215619300146

Moving Particle Semi-implicit method

Zidi Wang , ... Akifumi Yamaji , in Nuclear Power Plant Design and Analysis Codes, 2021

18.2.1 Governing equations

In the Lagrangian description of the Navier–Stokes equations, that is, mass conservation equation and momentum conservation equation, for an incompressible Newtonian fluid, can be written as follows:

(18.1) D ρ Dt = ρ · u = 0

(18.2) D u Dt = 1 ρ P + v 2 u + g + F ρ

where ρ is the density of the fluid, u is the flow velocity, P is the pressure, v is the kinematic viscosity, g is the gravitational acceleration, and F represents other forces, such as the surface tension.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128181904000188

CONTINUUM THEORY OF NONLINEAR VISCOELASTICITY

A.C. ERINGEN , R.A. GROT , in Mechanics and Chemistry of Solid Propellants, 1967

Lagrangian Description

In the Lagrangian description the surface of the body remains fixed. The stress and heat flux are respectively prescribed by:

(13.16) T K L N K x l , L = t ( n ) l C 1 K L N K N L j

(13.17) Q K N K = q ( n ) C 1 K L N K N L j

where N K is the exterior normal to the surface b(x) = 0. The factor on the right side of (13.16) and (13.17) is due to the area change caused by the deformation.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9781483198378500144

Passive tracer clustering and diffusion in random hydrodynamic flows

V.I. Klyatskin , in Dynamics of Stochastic Systems, 2005

Problem 66

Derive the Fokker–Planck equation for the gradient of particle concentration in divergence-free random velocity field in the case of the two-dimensional average linear shear u0 (r, t) = αy l, where r = (x, y), l = (1, 0).

Instruction

In the Lagrangian description, the particle concentration gradient in divergence-free velocity field satisfies the equations

d d t p x ( t | r 0 ) = p k ( t | r 0 ) u k ( r , k ) x , d d t p y ( t | r 0 ) = α p x ( t | r 0 ) p k ( t | r 0 ) u k ( r , k ) y , p ( 0 | r 0 ) = p 0 ( r 0 ) = ρ 0 ( r 0 ) .

Solution

(11.61) t P ( t , p | r 0 ) = { α p x p y + 1 8 D s ( 3 2 p 2 p 2 2 2 p k p l p k p l ) } P ( t , p | r 0 ) , P ( 0 , p | r 0 ) = δ ( p p 0 ( r 0 ) ) .

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780444517968500124

Introduction

Seiichi Koshizuka , ... Takuya Matsunaga , in Moving Particle Semi-implicit Method, 2018

Abstract

Particle methods are based on Lagrangian description and meshless discretization to simulate continuum mechanics. Advantages and disadvantages of the particle methods are discussed from these viewpoints. Moving particle semi-implicit (MPS) method is one of the particle methods. The fundamental idea of the MPS method is a weighted difference without the mesh. Particle interaction models are prepared for differential operators and substituted into the governing equations. Pressure Poisson equation is constructed to solve the pressure field implicitly, while the other terms are explicitly calculated. Differences between the MPS method and smoothed particle hydrodynamics (SPH) method are explained. The research history of the particle methods for continuum mechanics is summarized.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128127797000011

GENERALIZED GOVERNING EQUATIONS FOR MULTIPHASE SYSTEMS: AVERAGING FORMULATIONS

Amir Faghri , Yuwen Zhang , in Transport Phenomena in Multiphase Systems, 2006

4.2.2 Lagrangian Averaging

Lagrangian averaging is directly related to the Lagrangian description of a system, which requires tracking the motion of each individual fluid particle. Therefore, Lagrangian averaging is a very useful tool when the dynamics of individual particles are of interest. To obtain Lagrangian time averaging, it is necessary to follow a specific particle and observe its behavior for a certain time interval. Then, the behavior of this particle is averaged over the time interval.

For a generalized function Φ = Φ(X, Y, Z, t), X, Y, and Z are material coordinates moving with the particle, and X, Y, Z are functions of the spatial coordinates x, y, z, and time t, i.e.,

X = X ( x , y , z , t ) , Y = Y ( x , y , z , t ) , Z = Z ( x , y , z , t )

The most widely used Lagrangian averaging is time averaging, where the time average of the function Φ in time interval of Δt is

(4.16) Φ ¯ = 1 Δ t Δ t Φ ( X , Y , Z , t ) d t

Lagrangian time averaging is performed for a distinct particle moving in the field; therefore, X, Y, and Z in the time interval Δt are not fixed in space. This focus on specific particles moving in space and time distinguishes Lagrangian averaging from Eulerian time averaging, which treats a fixed point in space relative to the reference frame. An example from daily experience will serve to illustrate this difference. In order to monitor traffic on the highway, the speed of all cars passing a point can be measured and averaged over a certain time interval – a case of Eulerian averaging. To catch an individual speeder, the police must follow the vehicle of interest to measure its speed as it moves in space over a certain time interval – a case of Lagrangian averaging.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B978012370610250009X

Extended Algorithms

Seiichi Koshizuka , ... Takuya Matsunaga , in Moving Particle Semi-implicit Method, 2018

3.4 Arbitrary Lagrangian-Eulerian

The particle methods are usually based on the Lagrangian description, where the variables located at the particle move with the fluid. The convection terms in the Navier-Stokes equations are not necessary to evaluate. On the other hand the grid methods are mainly based on the Eulerian description, where the variables located at the nodes are fixed in space. The convection terms are necessary to evaluate. Here the particles in the particle methods and the nodes in the grid methods have the same meaning that the variables, velocity components, pressure, temperature, etc. are located. If particles or nodes are moved arbitrarily the description is called arbitrary Lagrangian-Eulerian (ALE). In the ALE description the convection terms are also necessary to evaluate.

In the ALE description, we can consider three particle positions in each time step as shown in Fig. 3.7. The position of particle i at time step k is denoted by r i k . The velocity is updated from u i k to u i L after the explicit and implicit steps of the MPS algorithm in this time step. The particle moves to

Figure 3.7. Movement of a particle in Arbitrary Lagrangian-Eulerian (ALE).

(3.99) r i L = r i k + Δ t u i L ,

in the fully Lagrangian description. All variables f are moved with the particle keeping the same values, i.e.,

(3.100) f k + 1 ( r i L ) = f k ( r i k ) .

The new-time position of particle i is then arbitrarily determined as r i k+1 with a velocity of the particle u i c

(3.101) r i k + 1 = r i k + Δ t u i c .

If the new position is the same as that of fully Lagrangian motion

(3.102) r i k + 1 = r i L ,

(3.103) u i c = u i L ,

This is the standard MPS method based on the Lagrangian description. If the new position is the same as the old position at time step k

(3.104) r i k + 1 = r i k ,

(3.105) u i c = 0 .

This is the fully Eulerian description. If the new position is arbitrarily given, we need to evaluate the change of variables as

(3.106) f k + 1 ( r i k + 1 ) = f k + 1 ( r i L Δ t u i a ) ,

(3.107) u i a = u i L u i c .

To calculate Eq. (3.106), we need to obtain the spatial distribution of f at time step k+1. This is equivalent to the calculation of the convection term.

In general the particle movement in one time step is smaller than that of the particle distance so as to satisfy the Courant condition. Thus the calculation of Eq. (3.106) is an interpolation of the variables at the particles.

Yoon et al. (1999a) proposed a method to calculate the interpolation equivalent to the convection in the ALE description. This method was named "meshless advection using flow-directional local-grid (MAFL)" and combined with MPS. First, a one-dimensional local grid through particle i is generated as shown in Fig. 3.8 in the flow direction of u i a . Temporary grid points are given; in Fig. 3.8, three additional grid points i+1, i−1, and i−2 are set with the spacing of Δr. Values at the temporary grid points are evaluated by weighted averages of their neighboring areas Ω i+1, Ω i-1, and Ω i-2, respectively. The neighboring area is limited by a circle of radius r e and lines vertical to the grid.

Figure 3.8. Flow-directional local grid in meshless advection using flow-directional local-grid (MAFL).

The new-time value at particle i is calculated, for instance, by the QUICK scheme (Leonard, 1979) as

(3.108) f i k + 1 = f i L Δ t Δ r | u i a | ( 3 8 f i + 1 + 3 8 f i 7 8 f i 1 + 1 8 f i 2 ) .

The calculation of the convection term with the flow velocity of u i a is equivalent to evaluate the interpolated value at the upwind position of r i L Δ t u i a on the flow-directional local grid. This concept is called "semi-Lagrangian" and widely used in the finite difference method.

In MPS-MAFL the temporary grid is generated for each particle at each time step, so that the global mesh is not necessary and mesh tangling does not occur. Since the temporary grid is one-dimensional, the calculation is not complicated even in three dimensions. Furthermore the knowledge of the finite difference schemes can be applied to the temporary grid. Thus MPS-MAFL is practical in the ALE description without losing the advantages of the particle methods.

In MPS-MAFL the source term of the pressure Poisson equation is represented by the divergence of the velocity because the particles are arbitrarily located and the particle number density is not proportional to the fluid density:

(3.109) u i = d n i j i [ ( u j u i ) ( r j r i ) | r j r i | 2 w ( | r j r i | ) ] ,

(3.110) 2 P k + 1 i = ρ Δ t u i ,

(3.111) u i L u i * Δ t = 1 ρ P k + 1 i .

A calculation example of flow-induced sloshing is shown in Fig. 3.9. Sloshing in a rectangular tank is excited by the flow entering from the left side and going out from the bottom (Yoon et al., 1999a). The inlet and outlet boundaries are calculated by the fixed particles (Eulerian), while the oscillating free surface is calculated by the moving particles (almost Lagrangian). The Eulerian treatment is simpler on the inlet and outlet boundaries than the Lagrangian treatment where the new particles are generated in the inlet and the outflowing particles are removed. However, the Lagrangian treatment is preferred on the moving free surface. The ALE treatment is particularly useful for this type of problems. The excitation of sloshing is affected by the water level and the inflow velocity. The calculated excitation condition showed good agreement with that of the experiment.

Figure 3.9. Calculation result of flow-induced sloshing in a rectangular tank using Moving Particle Semi-implicit-meshless advection using flow-directional local-grid (MPS-MAFL) (Yoon et al., 1999a).

One cycle of the bubble growth and departure in subcooled nucleate boiling was simulated using MPS-MAFL in two dimensions (Yoon et al., 2001). The result is shown in Fig. 3.10. A small bubble was initially located on the heated wall of which temperature was 110°C. The pressure was atmospheric and the bulk water temperature was 96°C. Steam in the bubble was assumed to be uniform and saturated. No particles were used in steam. The initial small bubble grew rapidly due to mass transfer from the thin temperature boundary layer near the hot wall on the bottom. The bubble detached from the bottom wall at around 0.03   s due to buoyancy force and rose up. Hot water rose up as well to follow the bubble motion. The bubble shrank due to condensation in the subcooled bulk water. The bubble growth and consequent heat transfer agreed well with the experiment.

Figure 3.10. Bubble growth and departure in nucleate boiling: heated wall temperature 110°C, saturation temperature 100°C, and bulk temperature 96°C (Yoon et al., 2001).

In this problem, small particles were placed in the thin boundary layers near the hot wall and the bubble surface, while larger particles were located in the bulk water where high spatial resolution is not necessary. The small particles stayed in the boundary layers using the ALE treatment. On the other hand it would not be possible in the fully Lagrangian treatment where the small particles move with their own velocities and some of them might flow out of the boundary layers. We can see that the number of particles increases when the bubble grows and rises up. More small particles are necessary for the boundary layer around the bubble. MPS-MAFL is applicable to both moving particles and new particles because the calculation procedure is the same for them. We need to interpolate a value at an arbitrary position from those at the existing particle positions. We do not need to care that the particle is moved or newly generated.

MPS-MAFL has been used in Yoon et al. (1999a, b, 2001), Koshizuka et al. (2000), and Heo et al. (2002).

Hu et al. (2017) proposed a new approach for the ALE particle method. The convection scheme is evaluated using LSMPS (Least Squares MPS) (Tamai et al., 2013; Tamai and Koshizuka, 2014) applied to the upwind half of the neighboring particles (Fig. 3.11). This can be expressed by the upwind weight function

Figure 3.11. Upwind weight function for Arbitrary Lagrangian-Eulerian (ALE) (Hu et al., 2017).

(3.112) w upwind ( r ) = { w ( | r i j | ) r i j n upwind , i > 0 0 otherwise ,

(3.113) r i j = r j r i ,

(3.114) n upwind , i = u i a / | u i a | ,

where r i and r j are position vectors of particles i and j, respectively. The proposed ALE particle method showed that accuracy was much improved than MPS-MAFL in a lid-driven cavity flow problem.

A particle shifting technique (Xu et al., 2009) was proposed for rearrangement of the particles to avoid clustering in Incompressible Smoothed Particle Hydrodynamics (ISPH) when the velocity divergence is employed as the source term in the pressure Poisson equation. After the particle shifting the variables at the particles are updated by using the first-order derivative and the distance vector, which is equivalent to evaluating the convection terms in the ALE approach.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128127797000035

Kinematics

In Fluid Mechanics (Fifth Edition), 2012

3.2 Particle and Field Descriptions of Fluid Motion

There are two ways to describe fluid motion. In the Lagrangian description, fluid particles are followed as they move through a flow field. In the Eulerian description, a flow field's characteristics are monitored at fixed locations or in stationary regions of space. In fluid mechanics, an understanding of both descriptions is necessary because the acceleration following a fluid particle is needed for the application of Newton's second law to fluid motion while observations, measurements, and simulations of fluid flows are commonly made at fixed locations or in stationary spatial regions with the fluid moving past the locations or through the regions of interest.

The Lagrangian description is based on the motion of fluid particles. It is the direct extension of single particle kinematics (e.g., see Meriam & Kraige, 2007) to a whole field of fluid particles that are labeled by their location, r o , at a reference time, t  = to . The subsequent position r of each fluid particle as a function of time, r(t;r o ,to ), specifies the flow field. Here, r o and t o are boundary or initial condition parameters that label fluid particles, and are not independent variables. Thus, the current velocity u and acceleration a of the fluid particle that was located at r o at time to are obtained from the first and second temporal derivatives of particle position r(t;r o ,to ):

(3.1) u = d r ( t ; r o , t o ) / d t and a = d 2 r ( t ; r o , t o ) / d t 2 .

These values for u and a are valid for the fluid particle as it moves along its trajectory through the flow field (Figure 3.4). In this particle-based Lagrangian description of fluid motion, fluid particle kinematics are identical to that in ordinary particle mechanics, and any scalar, vector, or tensor flow-field property F may depend on the path(s) followed of the relevant fluid particle(s) and time: F  = F[r(t;r o ,to ), t].

FIGURE 3.4. Lagrangian description of the motion of a fluid particle that started at location r o at time to . The particle path or particle trajectory r(t;r o ,to ) specifies the location of the fluid particle at later times.

The Eulerian description focuses on flow field properties at locations or in regions of interest, and involves four independent variables: the three spatial coordinates represented by the position vector x, and time t. Thus, in this field-based Eulerian description of fluid motion, a flow-field property F depends directly on x and t: F  = F(x, t). Even though this description complicates the calculation of a, because individual fluid particles are not followed, it is the favored description of fluid motion.

Kinematic relationships between the two descriptions can be determined by requiring equality of flow-field properties when r and x define the same point in space, both are resolved in the same coordinate system, and a common clock is used to determine the time t:

(3.2) F [ r ( t ; r o , t o ) , t ] = F ( x , t ) when x = r ( t ; r o , t o ).

Here the second equation specifies the trajectory followed by a fluid particle. This compatibility requirement forms the basis for determining and interpreting time derivatives in the Eulerian description of fluid motion. Applying a total time derivative to the first equation in (3.2) produces

(3.3) d d t F [ r ( t ; r o , t o ) , t ] = F r 1 d r 1 d t + F r 2 d r 2 d t + F r 3 d r 3 d t + F t = d d t F ( x , t ) when x = r ( t ; r o , t o ) ,

where the components of r are ri . In (3.3), the time derivatives of ri are the components ui of the fluid particle's velocity u from (3.1). In addition, ∂F/∂ri   = ∂F/∂xi when x  = r, so (3.3) becomes

(3.4) d d t F [ r ( t ; r o , t o ) , t ] = F x 1 u 1 + F x 2 u 2 + F x 3 u 3 + F t = ( F ) u + F t D D t F ( x , t ) ,

where the final equality defines D/Dt as the total time derivative in the Eulerian description of fluid motion. It is the equivalent of the total time derivative d/dt in the Lagrangian description and is known as the material derivative, substantial derivative, or particle derivative, where the final attribution emphasizes the fact that it provides time derivative information following a fluid particle.

The material derivative D/Dt defined in (3.4) is composed of unsteady and advective acceleration terms. (1) The unsteady part of DF/Dt, ∂F/∂t, is the local temporal rate of change of F at the location x. It is zero when F is independent of time. (2) The advective (or convective) part of DF/Dt, u⋅F, is the rate of change of F that occurs as fluid particles move from one location to another. It is zero where F is spatially uniform, the fluid is not moving, or u and ∇F are perpendicular. For clarity and consistency in this book, the movement of fluid particles from place to place is referred to as advection with the term convection being reserved for the special circumstance of heat transport by fluid movement. In vector and index notations, (3.4) is commonly rearranged slightly and written as

(3.5) D F D t F t + u F , or D F D t F t + u i F x i .

The scalar product u⋅F is the magnitude of u times the component of ∇F in the direction of u so (3.5) can then be written in scalar notation as

(3.6) D F D t F t + | u | s ,

where s is a path-length coordinate on the fluid particle trajectory x  = r(t;r o ,to ), that is, d r = e u d s with e u = u / | u | .

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780123821003100034

Response of composite panels to blast wave pressure loadings

K. Lee , S.W. Lee , in Blast Protection of Civil Infrastructures and Vehicles Using Composites, 2010

6.2.1 Non-linear solid shell element formulation for dynamic problems

The non-linear solid shell formulation is based on the total Lagrangian description that employs Green's strain and the second Piola–Kirchhoff stress. In the solid shell element formulation, the composite plate is treated as a three-dimensional solid, allowing thickness change and transverse shear deformation. The kinematics of deformation is described by vector variables only, not including any rotational angles. For dynamic problems, the mass matrix remains constant during the analysis. For geometrically non-linear static cases, it has been shown that large load increments can be used with the solid shell element formulation ( Park et al., 1995).

The assumed strain formulation is combined with the solid shell formulation to alleviate element locking. In this approach, an assumed strain field is carefully chosen independently of the displacement-dependent strain field. These two strain fields are related to each other by satisfying a set of compatibility equations. The assumed strain field within an element is expressed by assumed strain shape functions and their parameters that are eliminated at element level. Accordingly, these additional assumed strain parameters do not increase the number of unknowns at global level. Detailed description of the assumed strain solid shell formulation has been provided by Park et al. (1995) and Kim and Lee (1988).

For the dynamic formulation, the trapezoidal rule is chosen for the time integration. The trapezoidal rule requires iterations in order to obtain the dynamic equilibrium state at an instant in time. An incremental form of equilibrium and compatibility equations based on incremental expressions for stress and strain vectors are used. The blast wave is treated as dynamic pressure loading applied over the plate surface and the time variation of the blast pressure is approximated by the Friedlander decay function.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9781845693992500068

31st European Symposium on Computer Aided Process Engineering

M. Robles-Santacruz , ... A.R. Uribe-Ramírez , in Computer Aided Chemical Engineering, 2021

5 Conclusions.

In conclusion, the SPH method can properly and naturally represent the creation of bubbles due to its purely lagrangian description. The model reported in the methodology and solved with the SPH method generates representative results with good precision as is prove with the numerical results obtained in the validation test. The surface tension model allows a stable surface drop during the rise of the bubble; moreover, with this model it is possible to study the rise of drops for different regimes. According to the parameters used in simulations, the length of injections is the principal parameter that modify the size of bubbles. The size of bubbles don't be affected by the flows used in the simulations. The 2D simulation is enough to know the area and speed of the bubbles, however, it is not able to predict well the shape of the bubbles. So, a 3D simulation is preferable for cases where is necessary to know about the shapes, movement and other important events.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B978032388506550067X