A coupled finite volume and material point method for two-phase simulation of liquid-sediment and gas-sediment flows
AA coupled finite volume and material point method for two-phase simulationof liquid–sediment and gas–sediment flows
Aaron S. Baumgarten a , Benjamin L. S. Couchman a , Ken Kamrin b, ∗ a Department of Aeronautics and Astronautics, b Department of Mechanical Engineering,Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Abstract
Mixtures of fluids and granular sediments play an important role in many industrial, geotechnical, andaerospace engineering problems, from waste management and transportation (liquid–sediment mixtures) todust kick-up below helicopter rotors (gas–sediment mixtures). These mixed flows often involve bulk motion ofhundreds of billions of individual sediment particles and can contain both highly turbulent regions and static,non-flowing regions. This breadth of phenomena necessitates the use of continuum simulation methods, suchas the material point method (MPM), which can accurately capture these large deformations while alsotracking the Lagrangian features of the flow (e.g. the granular surface, elastic stress, etc.).Recent works using two-phase MPM frameworks to simulate these mixtures have shown substantialpromise; however, these approaches are hindered by the numerical limitations of MPM when simulatingpure fluids. In addition to the well-known particle ringing instability and difficulty defining inflow/outflowboundary conditions, MPM has a tendency to accumulate quadrature errors as materials deform, increasingthe rate of overall error growth as simulations progress. In this work, we present an improved, two-phasecontinuum simulation framework that uses the finite volume method (FVM) to solve the fluid phase equationsof motion and MPM to solve the solid phase equations of motion, substantially reducing the effect of theseerrors and providing better accuracy and stability for long-duration simulations of these mixtures.
Keywords: granular materials, poromechanics, FVM, MPM
1. Introduction
In this work, we are concerned with the numerical simulation of mixtures of fluids and grains in a widerange of applications and geometries. These mixtures are primarily composed of two materials: a porousgranular solid (e.g. sand, powder, soil, grains) and a fluid (e.g. water, gas, air). Mixtures of relativelydense fluids and sediments play an important role in many industrial and geotechnical engineering problems,from transporting large volumes of industrial wastes [1] to building earthen levees and dams [2]. On theother hand, mixtures of relatively light fluids (e.g. gases) and sand-like materials appear in vastly differentengineering problems, from understanding dust kick-up below helicopter rotors [3] to modeling locomotionof wheels on rough terrain [4].To solve these problems, engineers have traditionally relied on a myriad of empirical models developedover the last century (e.g. models for the effective viscosity of mixtures [5, 6, 7], the permeability of granularbeds [8], and the various regimes of slurry transport [1]). These empirical models are derived by couplingrelevant experimental observations to an understanding of the underlying physics governing the behaviorof these mixtures; however, such models can only describe specific regimes of mixed flow. To address an ∗ Corresponding author
Email address: [email protected] (Ken Kamrin)
Preprint submitted to Elsevier December 29, 2020 a r X i v : . [ c s . C E ] D ec ngineering problem that involves complex interactions of fluids and grains spanning many flow regimes orin a complex geometry, engineers require a more general modeling approach.A popular solution to this problem is modeling the microscopic behavior of these mixtures directly bysimulating each grain–grain collision and calculating the exact motion of the pore fluid numerically (e.g.LBM–DEM [9], CFD–DEM [10]). This type of approach is highly accurate and has been used to developcomplex constitutive relations for mixtures [11, 12, 13]; to predict wave generation from landslides [14]; andto model consolidation, shearing, and collapse of saturated soil columns [15, 16]. The primary limitationof this approach is difficulty simulating engineering problems involving large volumes of material. Sinceeach grain and its surrounding fluid must be tracked over time, simulations involving hundreds of billions ofsediment particles are computationally intractable. We therefore turn to a continuum modeling approach,where small scale phenomena are homogenized into bulk properties and behaviors.Recent work simulating fluid-grain mixtures as continua (see [17]) presents a versatile foundation, butsubstantial improvement is necessary before existing simulation frameworks can be used for the range ofengineering applications we seek to address. There are two primary challenges associated with such frame-works: i) accurate prediction of the stresses in the grains and fluid (i.e. the constitutive response of themixture; see [18, 19]) and ii) accurate numerical approximation of the equations of motion for mixed flows(see [20, 21]). It is this latter challenge that we seek to address in the present work.When represented as continua, the individual grains and fluid-filled pores of these mixtures are homog-enized into two independent material phases: the granular (or solid) phase and the fluid phase; each phasehas its own velocity ( v s and v f ), density ( ¯ ρ s and ¯ ρ f ), and internal state ( σ s , σ f , φ , etc.). This requires thattwo sets of mass, momentum, and energy conservation laws be satisfied simultaneously and that the motionof the phases be tracked independently (see [22, 23, 24, 25]). Further, these mixtures are known to undergoextremely large deformations (as in riverbed erosion; see [26]) with internal flows that are highly coupled tosurface topology changes (as in near-shore wave breaking; see [27]), requiring both Lagrangian-like featuretracking for the deforming surfaces and Eulerian-like evaluation of the governing equations to avoid the issuesassociated with entangling meshes on the interior.One continuum approach that has shown success in simulating these challenging problems is the two-phase material point method (MPM; see [20, 21, 28, 29]) which has been used to model water-jet–soilinteractions in [30], fluidization of granular beds in [31], submerged cone penetration in [32], slope failures in[19], and shear thickening of dense suspensions in [33]. A key advantage of MPM is that it naturally capturesLagrangian flow features by using persistent material point tracers while solving the equations of motion on atemporary Eulerian grid. Despite its proven robustness in modeling many solid-like materials through largedeformations, MPM is known to develop substantial integration errors when these deformations continue togrow, as occurs with pure fluids and mixtures with negligible solid volume fractions; these quadrature errorsincrease the rate of overall error growth during a simulation, limiting the use of MPM for long-duration fluidproblems (see [34, 35, 36, 37, 38]). Although substantial work has been done to reduce these errors in MPM(see [39, 40, 41, 42, 43, 44, 45]), it is desirable to avoid them completely where possible.The finite volume method (FVM) has also been applied to these types of problems, though usually insituations where both phases of the mixture are fully fluidized (e.g. sediment transport; see [46, 47, 48])or the solid phase remains static (e.g. fluid flow through subterranean faults and groundwater motion; see[49, 50]). A key advantage of FVM is that it naturally conserves mass, momentum, and energy by calculatingfluxes between neighboring cells and, since these cells do not move with the material, does not accumulatethe kind of quadrature errors observed in MPM. As a result, FVM can be applied to many fluid problemsthat are currently extremely difficult for MPM to model, such as turbulent flows in pipes and simulationsinvolving gases. FVM can also be augmented to track Lagrangian features of fluid flows numerically (see[51, 52, 53, 54]). Its long history of development and use for solving the equations of motion for fluids on bothEulerian and arbitrary Lagrangian-Eulerian grids (see [55, 56, 57, 58]) make it an ideal choice for replacingMPM when substantial material deformation is expected.In this work, we introduce a new approach for simulating fluid–sediment mixtures that combines theadvantages of MPM for simulating the granular phase of the mixture with the proven robustness of FVM forsimulating pure fluids. This approach discretizes the governing equations described in [21, 59] on two over-2apping simulation domains. Similar coupled approaches have shown success for fluid–structure interactionand solid–solid collision simulations (see [54, 60, 61, 62]). Here the granular phase is represented on a set ofmaterial point tracers with kinematic fields (e.g. velocity) defined on an associated finite element grid, andthe fluid phase is represented by average field quantities within a set of Eulerian finite volumes. In this way,we can use MPM to track the important Lagrangian features of the granular flow (e.g. solid volume fraction,elastic stress, granular free surface) and use FVM to solve the fluid equations of motion without sacrificingaccuracy in long-running turbulent flows.
2. Governing Equations
We begin describing our approach for simulating fluid–sediment mixtures by reviewing the system ofgoverning equations that we seek to solve. These equations derive from the continuum mixture theories of[22, 23, 24] (see Appendix A) and can also be found in [21] where they are simulated using the two-phaseMPM approach. In this formulation, we assume that the grains are rough (i.e. true contact can occurbetween grains; see [63]), incompressible with density ρ s , and essentially spherical with mean diameter d .Additionally, we assume that the grains are quasi-mono-disperse (no size-segregation during flow) and fullyimmersed in a compressible fluid having density ρ f and viscosity η .To capture the dynamic behavior of the coupled flows we are interested in, we consider the materials thatconstitute our mixture independently. The individual solid grains are homogenized into a single continuumbody called the solid (or granular ) phase , which has an effective density ¯ ρ s , velocity v s , and specific internalenergy ε s . Similarly, the fluid that fills the pore space between grains is homogenized into a single continuumbody called the fluid phase , which has an effective density ¯ ρ f , velocity v f , and specific internal energy ε f .If we let n denote the volume fraction of fluid in the mixture (also called the porosity of the mixture) and φ denote the volume fraction of solid (also called the granular packing fraction ) then the phase-wise effectivedensities ( ¯ ρ s and ¯ ρ f ) can be defined by the following relations, φ = 1 − n, ¯ ρ s = φρ s , ¯ ρ f = nρ f , (2.1)for all spatial points x in the mixed domain Ω .Since we are interested in tracking Lagrangian flow features during our simulations, it is useful to expressthe equations governing the behavior of the solid phase in their Lagrangian form. To do this, we define thesolid phase material derivative, D s /Dt , as follows, D s ψDt ≡ ∂ψ∂t + v s · ∇ ψ, (2.2)for an arbitrary scalar field ψ and with ∇ the spatial gradient operator (i.e. ∇ ≡ ∂/∂ x ).Mass conservation in the solid phase takes the following form, D s ¯ ρ s Dt = − ¯ ρ s div( v s ) , (2.3)with div() the spatial divergence operator (i.e. div( v s ) ≡ ∇ · v s ). Since we assume that ρ s is constant anduniform everywhere, this equation also allows us to track the evolution of the granular packing fraction andmixture porosity over time.Momentum conservation in the solid phase takes the following form, ¯ ρ s D s v s Dt = div( ˜ σ ) + ¯ ρ s g − f d − φ ∇ p f , (2.4)with ˜ σ the effective granular stress tensor, g the gravitational acceleration vector, f d the inter-phase dragvector, and p f the fluid phase pore pressure .The effective granular stress tensor, ˜ σ , has many different admissible forms depending both on the type ofmaterial used and the flow regimes under consideration. Commonly referenced examples include models for3et clays and dense sediments [64, 65], models for dry granular materials [66, 67], models for dense slurries[13, 68], and kinetic theories for fast moving granular gases [69, 70]. In this work, we assume that each ofthese models can be expressed generally as, ˜ σ = ˜ T (cid:0) ∇ v s , ¯ ξ (cid:1) , (2.5)with ˜ T a model-specific function of the solid phase velocity gradient and the current solid phase state vector ¯ ξ . This state vector fully describes the relevant properties of the granular material (e.g. deformation gradient,packing fraction, internal energy, fabric tensor, etc.) and can evolve over time. This equation of state oftentakes the place of an explicit energy conservation equation when isothermal or adiabatic assumptions aremade. In this work, we will primarily use the granular material models described in [21, 67].The inter-phase drag vector, f d , captures the tractions on particle surfaces that arise due to the relativemotion of the particles and pore fluid (e.g. Stokes drag, Darcy drag, etc.). This force (along with buoyancy)couples the behavior of the mixed phases. Following the works of [71, 72], we let the drag vector take theform, f d = 18 φ (1 − φ ) η d ˆ F ( φ, Re) ( v s − v f ) , (2.6)with Re ≡ ( nρ f d (cid:107) v s − v f (cid:107) ) /η (see [73]) and ˆ F ( φ, Re) a function given in those works.The fluid pore pressure, p f , is commonly given through an equation of state that has the following form, p f = ˆ p f (cid:0) ρ f , ε f (cid:1) , (2.7)with ˆ p f a model-specific function. In [20, 21] an isothermal, weakly-compressible, barotropic model is used;however, an ideal gas law or thermo-fluid models could also be chosen.Mass conservation in the fluid phase, expressed in an Eulerian frame, takes the following form, ∂ ¯ ρ f ∂t = − div(¯ ρ f v f ) . (2.8)Since we assume a compressible fluid fills the pores of the granular skeleton, we can use equations (2.1),(2.3), and (2.8) to uniquely determine the true fluid density, ρ f .Momentum conservation in the fluid phase, expressed in an Eulerian frame, takes the following form, ∂ ¯ ρ f v f ∂t = − div (cid:0) np f I + ¯ ρ f v f ⊗ v f (cid:1) + div( τ f ) + ¯ ρ f g + f d + p f ∇ n, (2.9)with τ f the effective fluid shear stress tensor and ⊗ the tensor product.The fluid shear stress tensor, τ f , is often expressed as, τ f = 2 η r D f , (2.10)with D f the deviatoric component of the fluid strain-rate tensor, D f ≡ ( ∇ v f + ∇ v (cid:62) f ) , and η r an effectivemixture viscosity. Although there is no universally accepted model for this component of stress across allfluid flow regimes (see [74]), it is generally agreed that the effective viscosity varies with the solid volumefraction φ ; first-order (see [5]), second-order (see [6]), and higher-order (see [24]) forms appear throughoutthe literature. As in [21], we will consider only the first-order form of this stress given in [5].Energy conservation in the fluid phase, expressed in an Eulerian frame, takes the following form, ∂ ¯ ρ f E f ∂t = − div (cid:0) (¯ ρ f E f + np f ) v f (cid:1) + div( τ f v f − q f ) + (¯ ρ f g ) · v f + (cid:0) p f ∇ n + f d (cid:1) · v s , (2.11)with q f the effective fluid heat flow vector and E f the specific total energy of the fluid with E f ≡ ε f + v f · v f .(Note that this form of the specific total energy does not include the contributions of microscopic fluid velocityfluctuations that may arise in mixed flows, such as from tortuosity; see [22, 75, 76, 77] and Appendix A forfurther discussion.)The effective fluid heat flow vector, q f , takes the following form, q f = − nk f ∇ ϑ f , (2.12)with k f the effective coefficient of thermal conductivity in the fluid and ϑ f the fluid phase temperature.4 . Mixed Weak-Form of Governing Equations In this section, we introduce a numerical approximation of the governing equations described above thatcan be evaluated using the two-phase, finite volume and material point method (FV-MPM). The governingequations are summarized as follows: D s ¯ ρ s Dt = − ¯ ρ s div( v s ) , ¯ ρ s D s v s Dt = div( ˜ σ ) + ¯ ρ s g − f d − φ ∇ p f ,∂ ¯ ρ f ∂t = − div(¯ ρ f v f ) ,∂ ¯ ρ f v f ∂t = − div (cid:0) ¯ ρ f v f ⊗ v f + np f I (cid:1) + div( τ f ) + ¯ ρ f g + f d + p f ∇ n,∂ ¯ ρ f E f ∂t = − div (cid:0) (¯ ρ f E f + np f ) v f (cid:1) + div( τ f v f − q f ) + (¯ ρ f g ) · v f + (cid:0) p f ∇ n + f d (cid:1) · v s . (3.1)Suppose we multiply the first two equations in (3.1) by a scalar test function, w , and a vector test function, w , respectively, and integrate both expressions over the same simulated domain, Ω . Suppose, further, thatwe integrate the final three equations in (3.1) over an arbitrary subdomain, Ω α , with boundary ∂ Ω α . Afterapplying the divergence theorem where relevant, the resulting set of equations can be separated into threedistinct groups: i) a weak expression of mass conservation in the solid phase to be evaluated on the materialpoints of MPM, (cid:90) Ω w D s ¯ ρ s Dt dv = − (cid:90) Ω w ¯ ρ s div( v s ) dv ; (3.2)ii) a weak expression of momentum conservation in the solid phase to be evaluated on the background finiteelement grid of MPM, (cid:90) Ω ¯ ρ s a s · w dv = (cid:90) Ω (cid:0) − ˜ σ : ∇ w + ¯ ρ s g · w − f d · w − p f I : ∇ ( n w ) (cid:1) dv + (cid:90) ∂ Ω τ s · w da ; (3.3)and iii) integral expressions for mass, momentum, and energy conservation in the fluid phase to be evaluatedon each subdomain of the finite volume grid of FVM, ∂∂t (cid:18) (cid:90) Ω α ¯ ρ f dv (cid:19) = − (cid:90) ∂ Ω α ¯ ρ f v f · ˆ n da,∂∂t (cid:18) (cid:90) Ω α ¯ ρ f v f dv (cid:19) = − (cid:90) ∂ Ω α (cid:0) ¯ ρ f v f ( v f · ˆ n ) + np f ˆ n (cid:1) da + (cid:90) ∂ Ω α τ f ˆ n da + (cid:90) Ω α (cid:0) ¯ ρ f g + f d + p f ∇ n (cid:1) dv,∂∂t (cid:18) (cid:90) Ω α ¯ ρ f E f dv (cid:19) = − (cid:90) ∂ Ω α (¯ ρ f E f + np f )( v f · ˆ n ) da + (cid:90) ∂ Ω α ( τ f v f − q f ) · ˆ n da + (cid:90) Ω α (¯ ρ f v f · g ) dv + (cid:90) Ω α (cid:0) p f ∇ n + f d (cid:1) · v s dv. (3.4)5Note that τ s in equation (3.3) denotes the prescribed traction on the solid phase at the simulated domainboundary ∂ Ω .)In order to evaluate these expressions numerically, we choose to represent the mixture fields on threeseparate solution spaces associated with the three groups above:i) a set of N m material point characteristic functions, { U p ( x , t ) | p ∈ [1 , N m ] } ,ii) a set of N n nodal basis functions, {N i ( x , t ) | i ∈ [1 , N n ] } ,iii) and a set of N v finite volume domains, { Ω α | α ∈ [1 , N v ] } ,with x a spatial position within the simulated domain Ω and t a point in time within the simulated period [0 , t end ] . A visual representation of these different spaces is shown in Figure 1 for general problems and inFigure 2 for the sand–air mixture simulation shown in Figure 8a. Figure 1: The basic method of solving mixture problems using the finite volume and material point method (FV-MPM). a)Define the mixture and initial configuration including densities, porosities, stresses. b) Define the initial solid and fluid phasecontinuum bodies ( B s and B f , respectively). c) Break the continuum bodies into piecewise-defined chunks of material: usematerial points to track the solid phase and finite volumes to track the fluid ( P p is the body part associated with the p thmaterial point). d) Define finite element nodal basis functions and simulation domain. e) Solve the discrete equations of motionfor the mixture using FV-MPM. Consider first the material points. Suppose the initial domain of the solid phase continuum, B s , iscomposed of N m distinct subdomains, {P p | p ∈ [1 , N m ] } , as in Figure 1c). We can then use these subdomainsto define the material points characteristic functions as follows, N m (cid:88) p =1 U p ( x ,
0) = (cid:26) x ∈ B s , else and U p ( x ,
0) = (cid:26) x ∈ P p , else ∀ p ∈ [1 , N m ] . (3.5)6 ) 𝑏) Figure 2: An example of a two-phase finite volume and material point method (FV-MPM) simulation highlighting the differentsolution spaces associated with each phase of the mixture. a) A snapshot from the simulation of a sand–air mixture with windblowing from left to right; the granular phase is represented by black material tracers and the fluid phase is colored according tothe local effective temperature, ϑ f . b) A region of this simulation showing the material point characteristic functions, U p ( x , t ) ,visualized using their GIMP representations (squares; see [78]); the nodal basis functions, N i ( x ) , and associated nodes (blackcircles); and the finite volume subdomains, Ω α (triangles). Here, we chose nodal basis functions and finite volume subdomainsthat are defined on the same triangular grid. During a simulation, the solid phase domain and its p th subdomain will evolve to B ts and P tp , respectively,and this evolution can be tracked using the characteristic functions if we require that they be co-moving withthe material (i.e D s ( U p ( x , t )) /Dt = 0 ). In practice, it is also useful to define a set of material point weights, v p ( t ) , and centroids, x p ( t ) , as follows, v p ( t ) = (cid:90) Ω U p ( x , t ) dv, x p ( t ) = 1 v p ( t ) (cid:90) Ω x U p ( x , t ) dv. (3.6)As is common in MPM, we let these characteristic functions define the solid phase material fields ( ¯ ρ s , ˜ σ ,and ¯ ξ ) as, ¯ ρ s = N m (cid:88) p =1 ¯ ρ sp ( t ) U p ( x , t ) , ˜ σ = N m (cid:88) p =1 ˜ σ p ( t ) U p ( x , t ) , ¯ ξ = N m (cid:88) p =1 ¯ ξ p ( t ) U p ( x , t ) , (3.7)with the test function w defined similarly, w = N m (cid:88) i =1 w p U p ( x , t ) . (3.8)Here ¯ ρ sp ( t ) , ˜ σ p ( t ) , ¯ ξ p ( t ) , and w p are the coefficients associated with the p th characterstic function and, byextension, the value of these fields at the p th material point centroid, x p ( t ) . Consider now the finite element nodal basis functions. This set of N n continuous ( C ) functions isdefined on an Eulerian computational grid, see Figure 1d, with the property that (cid:80) N n i =1 N i ( x ) = 1 . Since7hese functions are not co-moving with the material, we have, ∂ (cid:0) N i ( x ) (cid:1) /∂t = 0 . As in MPM, we use thesefunctions to define the solid phase kinematic fields, a s and v s , as follows, D s v s Dt = a s = N n (cid:88) i =1 a si ( t ) N i ( x ) , v s = N n (cid:88) i =1 v si ( t ) N i ( x ) , (3.9)along with the test function w and the mixture porosity field n , w = N n (cid:88) i =1 w i N i ( x ) , n = N n (cid:88) i =1 n i ( t ) N i ( x ) . (3.10)Here a si ( t ) , v si ( t ) , w i , and n i ( t ) are the coefficients associated with the i th nodal basis function but notnecessarily the value of these fields at the position i th node (e.g. B-spline coefficients; see [79]). The final space on which our mixture solutions are defined is the set of N v non-overlapping finite volumes.This set of volumes is defined on an Eulerian computational grid, which may be unique from the grid associatedwith the N n nodal bases , with ∪ N v α =1 Ω α = Ω and Ω α the domain encompassed by the α th finite volume, seeFigure 1c. Each of these volumes has an associated boundary ∂ Ω α , an associated weight, V α , and anassociated centroid, X α , defined as, V α ≡ (cid:90) Ω α dv, X α ≡ (cid:90) Ω α x dv, (3.11)which are fixed in the Eulerian frame, such that ∂ (cid:0) V α (cid:1) /∂t = 0 and ∂ (cid:0) X α (cid:1) /∂t = . As is common in FVM,we track the fluid phase effective density ( ¯ ρ f ), momentum ( ¯ ρ f v f ), and energy ( ¯ ρ f E f ) by defining cell-wiseaverages for each of these fields over each of the N v finite volumes: (cid:104) ¯ ρ f (cid:105) α ( t ) = 1 V α (cid:90) Ω α ¯ ρ f dv, (cid:104) ¯ ρ f v f (cid:105) α ( t ) = 1 V α (cid:90) Ω α ¯ ρ f v f dv, (cid:104) ¯ ρ f E f (cid:105) α ( t ) = 1 V α (cid:90) Ω α ¯ ρ f E f dv, (3.12)with (cid:104) ¯ ρ f (cid:105) α ( t ) , (cid:104) ¯ ρ f v f (cid:105) α ( t ) , and (cid:104) ¯ ρ f E f (cid:105) α ( t ) are time-dependent cell averages associated with the volume Ω α . We now return to the weak- and integral-form governing equations defined in (3.2), (3.3), and (3.4).Substitution of the solution fields described above into these equations can be simplified by using the mappingmatrices [ S ] , [ G ] , [ A ] , [ B ] , and [ M ] with components defined as follows, S ip = 1 v p (cid:90) Ω N i ( x ) U p ( x , t ) dv, G ip = 1 v p (cid:90) Ω ∇N i ( x ) U p ( x , t ) dv, A iα = 1 V α (cid:90) Ω α N i ( x ) dv, B ij = (cid:90) Ω N i ( x ) N j ( x ) dv, M ij = N m (cid:88) p =1 (cid:90) Ω ¯ ρ sp U p ( x , t ) N i ( x ) N j ( x ) dv. (3.13) [ S ] and [ G ] are the standard, time-varying MPM mapping matrices that allow for cross integration of fieldsbetween the material points and the finite element nodes, [ A ] is the FE-FVM mapping matrix that allows8or integration of fields defined on the finite element bases over the set of finite volumes, [ B ] is the volumematrix of the finite element grid, and [ M ] is the time-dependent finite element mass matrix associated withthe solid phase momentum balance equation.Using these matrices and equations (3.5), (3.7), (3.8), and (3.9), we can re-express the solid phase massconservation expression in (3.2) as follows, ddt (cid:0) ¯ ρ sp (cid:1) = − N n (cid:88) i =1 ¯ ρ sp v si · G ip , ∀ p ∈ [1 , N m ] , (3.14)though it is sometimes more useful to express this equation in terms of the fixed material point masses, m p = v p ¯ ρ sp , with, ddt (cid:0) m p ) = 0 , v p ddt (cid:0) v p (cid:1) = N n (cid:88) i =1 v si · G ip , ∀ p ∈ [1 , N m ] . (3.15)Similarly, recognizing that the components of w can be defined arbitrarily, we can substitute (3.7), (3.9),and (3.10) into the solid phase momentum conservation expression in (3.3) to find, N n (cid:88) j =1 M ij a sj = f int i + f ext i + f drag i + f buoy i + f τ i , ∀ i ∈ [1 , N n ] f int i = − N m (cid:88) p =1 G ip (cid:62) v p ˜ σ p , f ext i = N m (cid:88) p =1 S ip m p g , f drag i = − (cid:90) Ω f d N i ( x ) dv, f buoy i = N n (cid:88) j =1 (1 − n j ) (cid:90) Ω ∇ (cid:0) N i ( x ) N j ( x ) (cid:1) p f dv, f τ i = (cid:90) ∂ Ω N i ( x ) τ s da. (3.16) f int i , f ext i , and f τ i are the common MPM nodal force vectors and can be calculated efficiently once thetime-varying coefficients of the [ S ] and [ G ] matrices are determined. f drag i and f buoy i are new coupling termsbetween the MPM and FVM domains and, in practice, require more careful calculation. In particular, thefluid pore pressure p f and fluid phase velocity v f , which are necessary to calculate these terms, must bereconstructed from the N v finite volume cell average fluid properties according to the procedure describedin the next section. Some useful approximations of these terms can also be found in Appendix B.Finally, consider the integral-form system of equations described in (3.4). It is useful, particularly forevaluating the fluxes across volume boundaries, to re-express this system of equations using the fluid state9ector, u = { ρ f ; ρ f v f ; ρ f E f } , as follows, ddt (cid:32) (cid:104) ¯ ρ f (cid:105) α (cid:104) ¯ ρ f v f (cid:105) α (cid:104) ¯ ρ f E f (cid:105) α (cid:33) = F int α + F ext α + F drag α + F buoy α , ∀ α ∈ [1 , N v ] , F int α = − V α (cid:90) ∂ Ω α ˆ f ( u , n ; ˆ n ) da, F ext α = (cid:32) (cid:104) ¯ ρ f (cid:105) α g (cid:104) ¯ ρ f v f (cid:105) α · g (cid:33) , F drag α = 1 V α (cid:90) Ω α N n (cid:88) i =1 (cid:32) f d f d · v si (cid:33) N i ( x ) dv, F buoy α = 1 V α (cid:90) Ω α N n (cid:88) i =1 N n (cid:88) j =1 (cid:32) − n i ) p f (1 − n i ) p f v sj (cid:33) ∇N i ( x ) N j ( x ) dv. (3.17)with ˆ f the oriented flux vector associated with the outward pointing face normal, ˆ n , defined as follows, ˆ f ( u , n ; ˆ n ) = n (cid:32) ρ f v f · ˆ n ρ f v f ( v f · ˆ n ) + p f ˆ n ( ρ f E f + p f )( v f · ˆ n ) (cid:33) − (cid:32) τ f ˆ n ( τ f v f − q f ) · ˆ n (cid:33) . (3.18)Numerical approximation of the finite volume forces in (3.17) is possible using numerical quadrature, but itis also necessary to reconstruct the fluid phase fields within each cell. Since this reconstruction may not becontinuous across the cell boundaries, ∂ Ω α , we introduce a numerical flux function, ˆ h ( u + , u − , n ; ˆ n ) , as in[57], such that, F int α ≈ − V α (cid:90) ∂ Ω α ˆ h ( u + , u − , n ; ˆ n ) da. (3.19)Here u + and u − represent the two reconstructed fluid state vector values on each side of the boundary ∂ Ω α .In this work, we implement the entropy adjusted Roe flux function (see [55, 80]), which has the followingform, ˆ h ( u + , u − , n ; ˆ n ) = (cid:2) ˆ f ( u + , n, ˆ n ) + ˆ f ( u − , n, ˆ n ) − | A | ( u + − u − ) (cid:3) , (3.20)where | A | is the positive definite matrix formed from the flux Jacobian, ∂ ˆ f /∂ u . To implement the entropyfix, we substitute ¯ A for A as follows, ( u + − u − ) = m (cid:88) k =1 b k r k , (cid:12)(cid:12) ¯ A (cid:12)(cid:12) ( u + − u − ) = m (cid:88) k =1 ˆ Q ( λ k ) b k r k , (3.21)with { λ k } and { r k } the sets of eigenvalues and right eigenvectors of A , { b k } defined by the relationshipabove, and ˆ Q ( λ k ) defined as: ˆ Q ( λ k ) = (cid:40) | λ k | if | λ k | ≥ δ a , (cid:0) ( λ k /δ a ) + δ a ) if | λ k | < δ a , (3.22)with δ a = 0 . a for eigenvalues associated with acoustic waves and δ a = 0 for eigenvalues associated withadvective waves ( ¯ a is the Roe averaged speed of sound). (Note that | ¯ A | and | A | are generally identical forthe cases considered in this work with the exception of section 8.6.)In this way, we can fully describe the mixture fields and their evolution in terms of the discrete coefficients, ¯ ρ sp ( t ) , ˜ σ p ( t ) , ¯ ξ p ( t ) , a si ( t ) , v si ( t ) , n i ( t ) , (cid:104) ¯ ρ f (cid:105) α ( t ) , (cid:104) ¯ ρ f v f (cid:105) α ( t ) , (cid:104) ¯ ρ f E f (cid:105) α ( t ) , (3.23)and the time-dependent deformations of the material point characteristic functions (further discussion below).10 . Reconstructing Fluid Fields within Finite Volumes As noted above, in order to evaluate the integral expressions in (3.16) and (3.17), we need to approximatethe values of the fluid fields ( p f , ¯ ρ f , ρ f , v f , E f , τ f , and q f ) at all points in the simulated domain usingthe set of N v finite volume cell averages and the N n porosity coefficients. To do this, we reconstruct the true fluid phase density ( ρ f ), momentum ( ρ f v f ), and energy ( ρ f E f ) fields using the second-order approachdescribed in [57].We begin by defining the average porosity of each finite volume, (cid:104) n (cid:105) α , as follows, (cid:104) n (cid:105) α = N n (cid:88) i =1 A iα n i , ∀ α ∈ [1 , N v ] (4.1)and approximate the cell-wise averages of these true fluid fields as, (cid:104) ρ f (cid:105) α ≈ (cid:104) ¯ ρ f (cid:105) α (cid:104) n (cid:105) α , (cid:104) ρ f v f (cid:105) α ≈ (cid:104) ¯ ρ f v f (cid:105) α (cid:104) n (cid:105) α , (cid:104) ρ f E f (cid:105) α ≈ (cid:104) ¯ ρ f E f (cid:105) α (cid:104) n (cid:105) α , ∀ α ∈ [1 , N v ] . (4.2)If we take these cell-wise averages to be the value of these fields at the cell centroids, then we can approximatethe fluid phase fields within each volume as follows, ρ f ≈ (cid:104) ρ f (cid:105) α + (cid:104)∇ ρ f (cid:105) α · ( x − X α ) ∀ x ∈ Ω α ,ρ f v f ≈ (cid:104) ρ f v f (cid:105) α + (cid:104)∇ ρ f v f (cid:105) α ( x − X α ) ∀ x ∈ Ω α ,ρ f E f ≈ (cid:104) ρ f E f (cid:105) α + (cid:104)∇ ρ f E f (cid:105) α · ( x − X α ) ∀ x ∈ Ω α , (4.3)with (cid:104)∇ ρ f (cid:105) α , (cid:104)∇ ρ f v f (cid:105) α , and (cid:104)∇ ρ f E f (cid:105) α the reconstructed field gradients in Ω α . To estimate these gradients,we use the flux limited reconstruction described by Barth & Jespersen in [57], (cid:104)∇ ρ f (cid:105) α = Φ ρ f [ ∇ ρ f ] α , (cid:104)∇ ρ f v f (cid:105) α = Φ v f [ ∇ ρ f v f ] α , (cid:104)∇ ρ f E f (cid:105) α = Φ E f [ ∇ ρ f E f ] α , (4.4)with [ ∇ ρ f ] α , [ ∇ ρ f v f ] α , and [ ∇ ρ f E f ] α the gradients associated with a linear, least squares fit to the sur-rounding centroid data and Φ ρ f , Φ v f , and Φ E f the maximum admissible values in the range [0 , that ensurethe local field reconstructions are bounded by the surrounding cell maxima/minima. With the true fluidfields and their gradients determined, we can approximate the effective fluid density field as, ¯ ρ f ≈ ρ f N n (cid:88) i =1 n i N i ( x ) , (4.5)and use appropriate equations of state (e.g. the ideal gas law or a barotropic fluid model; see [21, 55]) tocalculate the remaining fields ( p f , τ f , and q f ). Further discussion of this approach, along with potentialimprovements that explicitly account for local variations in the porosity field, can be found in Appendix B.
5. Convection of Material Points
One of the primary advantages of MPM is the ability to track Lagrangian flow features and history-dependent material properties on the material point characteristic functions. Since the material pointsthemselves move with the material in a Lagrangian frame (see Figure 3), the convection of these additionalproperties occurs naturally and does not need to be explicitly calculated (i.e. there is no need to introduceadditional advection equations). However, the time-dependent nature of the characteristic functions requiresthat the mapping matrices [ S ] and [ G ] be recalculated frequently and that this motion be numericallyapproximated, both of which come at significant computational cost. In this work, we approximate the11otions and deformations of the material points using the original MPM approach from [28] in section 8.3and the uGIMP approach from [39] in sections 8.1, 8.2, 8.4, 8.5, and 8.6. (uGIMP provides better stabilitywhen linear basis functions are chosen.) Both of these approaches track the changes to the material pointweights as in equation (3.15) and centroids as follows, ddt (cid:0) x p (cid:1) ≈ N n (cid:88) i =1 S ip v si ∀ p ∈ [1 , N m ] . (5.1)Although higher-order approximations exist in the literature (e.g. GIMP [78], CPDI [39], and CPDI2 [41]),these methods are known to be unstable in problems involving extremely large deformations. Figure 3: Schematic showing a representative finite element grid and material point discretization. The continuum body andthe p th material subdomain are shown in their a) initial configurations ( B and P p ) and b) configurations at time t ( B t and P tp ). The p th material point characteristic function, U p ( x , t ) , is determined by the position of the p th material subdomain. With the material point motion determined, we eliminate the need for an additional convection equationfor the porosity field, n , by defining a weak, mass conserving mapping between the material point densityfield and finite element porosity field using the test function w as follows, (cid:90) Ω (1 − n ) ρ s w dv = (cid:90) Ω ¯ ρ s w dv, (5.2)which is equivalently, N n (cid:88) j =1 B ij (1 − n j ) ρ s = N m (cid:88) p =1 S ip m p , ∀ i ∈ [1 , N n ] . (5.3)Similarly, to convect solid phase momentum, we introduce an approximation of the velocity field, v ∗ s ( x , t ) ,that is defined on the material point characteristic functions, v ∗ s ( x , t ) = N p (cid:88) p =1 v ∗ sp ( t ) U p ( x , t ) , (5.4)and can be tracked independently from the true velocity field, v s ( x , t ) . We then define a weak, momentumconserving mapping between these fields using the test function w as follows, (cid:90) Ω ¯ ρ s v s · w dv = (cid:90) Ω ¯ ρ s v ∗ s · w dv, (5.5)which is equivalently, N n (cid:88) j =1 M ij v sj = N p (cid:88) p =1 S ip m p v ∗ sp , ∀ i ∈ [1 , N n ] . (5.6)12n this way, we can define the coefficients associated with the node-based solid velocity field and mixtureporosity field at any time using the persistent material point densities and velocities. The evolution of thematerial point densities is defined in (3.15) while the evolution of the material point velocities is defined as, ddt (cid:0) v ∗ sp (cid:1) = N n (cid:88) i =1 S ip a si . (5.7)This approach is sometimes referred to as the FLIP projection method (see [81, 82]), though there are severalother valid projections that appear in the MPM literature (e.g PIC [28], APIC [83], and XPIC [84]).
6. Numerical Algorithm
In order to find time-accurate solutions to the system of equations in (3.1), we need to describe thetime-history of the solution coefficients in (3.23) in some consistent manner. Toward this end, we discretizethe time domain, [0 , t end ] , into a set of N s + 1 discrete time markers, { t k ∀ k ∈ [0 , N s ] } , with, t k = k ∆ t and ∆ t = t end /N s . (6.1)We can then denote the values of our solution coefficients at these time markers using the convention, ψ kj = ψ j ( t k ) , with ψ j ( t ) an arbitrary, time-dependent function, and determine their values using the numericalalgorithm described below. As in [28], calculations in this numerical algorithm can be simplified by usingdiagonalized computational matrices. Here [ M ] k and [ B ] are diagonalized by summing each row into [ M D ] k and [ B D ] . This substitution introduces some error but still produces a consistent method (see [85]).Suppose the following parameters are known at t k : (i) the mapping matrices, [ S ] k , [ G ] k , [ A ] , [ B D ] , and [ M D ] k ; (ii) the solution coefficients, ¯ ρ ksp , ˜ σ kp , ¯ ξ kp , (cid:104) ¯ ρ f (cid:105) kα , (cid:104) ¯ ρ f v f (cid:105) kα , and (cid:104) ¯ ρ f E f (cid:105) kα ; (iii) the coefficients of thematerial point approximation of the velocity field, v ∗ ksp ; and (iv) the numerical representation of the materialpoint characteristic functions, U p ( x , t k ) , (e.g. v kp and x kp for basic MPM [28] and uGIMP [39]). We canuse a first-order, Forward Euler time integration method (with the update-stress-last MPM approach from[35, 86]) to determine the value of the solution coefficient at t k +1 according to the following steps:(1) Determine the mixture porosity coefficients, as in (5.3), using the diagonal matrix [ B D ] : B Dii (1 − n ki ) ρ s = N m (cid:88) p =1 S kip m p ∀ i ∈ [1 , N n ] . (6.2)(2) Determine the solid phase velocity coefficients, as in (5.6), using the diagonal matrix [ M D ] : M Dkii v ksi = N p (cid:88) p =1 S kip m p v ∗ ksp , ∀ i ∈ [1 , N n ] . (6.3)(3) Calculate the nodal force vectors, ( f int i ) k and ( f ext i ) k , as in (3.16).(4) Construct an approximation of the fluid phase fields, ρ f , v f , E f , and ¯ ρ f , within each finite volumecell, as in (4.4).(5) Use the reconstructed fluid phase fields to calculate the fluid phase pressures and stresses, p f and τ f ; the fluid phase heat flux, q f ; the nodal force vectors, ( f drag i ) k and ( f buoy i ) k , as in (3.16); andthe finite volume force vectors, ( F int α ) k , ( F ext α ) k , ( F drag α ) k , and ( F buoy α ) k , as in (3.17) and (3.19).Numerical quadrature can be particularly helpful in these calculations. (See Appendix B for alternativeapproximations of the drag and buoyancy terms.)136) Use an explicit time integrator to obtain the updated fluid phase field coefficients from the equationsof motion: (cid:104) ¯ ρ f (cid:105) k +1 α (cid:104) ¯ ρ f v f (cid:105) k +1 α (cid:104) ¯ ρ f E f (cid:105) k +1 α = (cid:104) ¯ ρ f (cid:105) kα (cid:104) ¯ ρ f v f (cid:105) kα (cid:104) ¯ ρ f E f (cid:105) kα + ∆ t (cid:2) ( F int α ) k + ( F ext α ) k + ( F buoy α ) k + ( F drag α ) k (cid:3) , ∀ α ∈ [1 , N v ] . (6.4)(7) Determine the solid phase acceleration coefficients associated with the finite element grid nodes fromthe equations of motion and appropriate boundary conditions: M Dkii a ksi = ( f int i ) k + ( f ext i ) k + ( f drag i ) k + ( f buoy i ) k + f τ i , ∀ i ∈ [1 , N n ] . (6.5)(8) Use an explicit time integrator to obtain the updated solid phase velocity coefficients associated withthe finite element grid nodes; treat the grid as though it were Lagrangian: v k +1 (cid:48) si = v ksi + ∆ t a ksi , ∀ i ∈ [1 , N n ] . (6.6)(9) Update the solid phase material state vector, ¯ ξ p , and effective granular stress, ˜ σ p , according to therelevant constitutive update procedure. For stability, L p , the average material point velocity gradientis often used: L k +1 p = N n (cid:88) i =1 v k +1 (cid:48) si ⊗ G kip , ˜ σ k +1 p = T (cid:0) L k +1 p , ¯ ξ k +1 p (cid:1) , ∀ p ∈ [1 , N m ] . (6.7)(10) Map the solid phase nodal accelerations and velocities to the material points to update their velocityapproximations, positions, and densities: v ∗ k +1 sp = v ∗ ksp + ∆ t N n (cid:88) i =1 S kip a ksi , ∀ p ∈ [1 , N m ] , x k +1 p = x kp + ∆ t N n (cid:88) i =1 S kip v k +1 (cid:48) si , ∀ p ∈ [1 , N m ] ,v k +1 p = v kp (cid:18) t N n (cid:88) i =1 G kip · v k +1 (cid:48) si (cid:19) , ∀ p ∈ [1 , N m ] , ¯ ρ k +1 sp = m p /v k +1 p , ∀ p ∈ [1 , N m ] . (6.8)(11) Update the diagonal mass matrix, [ M D ] k +1 , and mapping matrices, [ S ] k +1 and [ G ] k +1 , according totheir definitions in (3.13). Numerical quadrature can be particularly helpful in these calculations,but care must be taken in updating the location and weights of quadrature points at this step (see[39, 41, 78] for further discussion).(12) Increment k to k + 1 ; go to (1).This algorithm is consistent with the governing equations in (3.1) and is used to generate the results insections 8.1 and 8.2. An alternative numerical algorithm is described in Appendix E and is used to generatethe results in sections 8.3, 8.4, 8.5, and 8.6. This alternative algorithm is more complex but provides betteraccuracy and stability of the simulated fluid. Both numerical algorithms have two important conditionsthat must be satisfied to ensure their stability. First, the Forward Euler time-integration method hasexplicit expressions for each coefficient update; this requires consideration of the Courant–Friedrichs–Lewy(CFL; see [87]) condition for numerical methods (i.e. the time increment ∆ t must be small enough that thenumerical stencil for any explicit equation captures the motion of physical waves in the real system during14hat increment). In practice, we find the most limiting wave-speed to be the acoustic waves in the solid orfluid domains (depending on the specific elastic moduli and material properties of the problem). Second,the numerical evaluation of the fluid phase shear stresses and heat fluxes at the finite volume boundaries aswell as the calculation of the inter-phase drag terms, ( f drag i ) k and ( F drag α ) k , can be unstable if the chosentime increment is not small enough to capture the decay rates associated with each of these terms; theseinstabilities are not unique to this framework and are common to many computational codes for transientproblems. Full discussion of these stability conditions, an analysis of the overall accuracy of this simulationapproach, and alternative numerical procedures can be found in Appendices C, D, and E.
7. Order of Accuracy
Suppose the numerical algorithm described above is implemented for a particular problem using timeincrement ∆ t , a finite volume mesh with characteristic length h α (e.g. the average distance between cellcenters), a finite element grid with characteristic length h i (e.g. the average distance between grid nodes),and a set of material points with characteristic length h p (e.g. the average distance between the initialpositions of the material point centroids). If the stability criteria mentioned above are satisfied and equations(3.14), (3.16), (3.17), and (5.3) were evaluated exactly, the results of the algorithm would have numericalerrors that were proportional to the truncation errors in the time-integration method (e.g. Forward Euler, O (∆ t ) ) and the discretization errors determined by the choices of finite volume mesh, nodal basis functions,and material point characteristic functions (e.g. Barth & Jesperson reconstructions, O ( h α ) ; non-overlappingmaterial point characteristics, O ( h p ) ; l th-order nodal bases, O ( h li ) ). In practice, however, the use of commonMPM quadrature schemes (e.g. uGIMP, CPDI, etc.), diagonal mapping matrices [ M D ] and [ B D ] , and theapproximations of the inter-phase forces described in Appendix B, limits the order of convergence of theoverall method to the order of the truncation and quadrature errors associated with each of these terms: O ( h p ) + O ( h p /h i ) + O ( h i ) + O ( h α ) + O (∆ t ) . The O ( h p /h i ) term is common to particle-in-cell (PIC)methods and is usually reduced by increasing the number of material points per grid cell . (Note that inthe special case of a regular Cartesian grid with B-spline basis functions, the convergence rate becomes O ( h p ) + O ( h p /h i ) + O ( h i ) + O ( h α ) + O (∆ t ) ; see [79].) A full analysis of these convergence rates can befound in Appendix C.
8. Numerical Examples
Models for flow through porous media have practical importance for a number of engineering applications,including predicting groundwater movement and gas seepage through packed sediments. Laws for the motionof fluids through porous solids have been developed and generalized for more than a century (see [8, 71, 72,73, 88]) and have been implemented in many numerical frameworks (e.g. [49]). In the case of uniaxial flowor for flows along a streamline, the apparent fluid velocity u (i.e. the fluid volume flux divided by the totalflow area, u = q/A ) can be expressed as a function of the soil permeability K , fluid pressure change ∆ p f ,and the length of the soil section L : u = K ∆ p f L . (8.1)Here K is an empirical measure of resistance due to drag on the fluid as it flows through the pores of thesolid.We test our numerical framework on an example problem of this type. Consider a 2m long, 10cm tallpipe filled with water ( κ = 2 . GPa, ρ f = 1000 kg/m , and η = 1 mPa · s) and a L = 1 m section of packedgranular solid as in Figure 4a. The porous solid in this problem can be modeled as a linear elastic materialwith effective Young’s modulus E = 10 MPa, effective Poisson’s ratio ν = 0 . , true density ρ s = 2650 kg/m ,solid volume fraction φ ∈ [0 . , . , and average grain diameter d = 1 mm. (Note that the porous solid canbe compressed without violating the assumption that the grains themselves are incompressible.) If we modelthe drag function from equation (2.6) using the Carman–Kozeny formula, F ( φ, Re) = 10 φ/ (1 − φ ) , then15he permeability of this solid section can be expressed as K = d (1 − φ ) / η φ . With this determined,we can express the true fluid velocity in the pipe, v fx ( x ) , analytically in terms of the volume fraction of thesolid, φ , and the local porosity in the pipe, n ( x ) : v fx ( x ) = 1 n ( x ) d η (1 − φ ) φ ∆ p f L . (8.2)Given this exact solution to the governing equations in (3.1), our numerical framework can be evaluated bycomparing its predicted solutions with those from (8.2).
Figure 4: Numerical test of flow through a porous solid: a) schematic diagram of test; b) finite volume (FV) grid, finite element(FE) grid, and material points (MP) used in numerical test; c) a representative analytical (dashed line) and simulated (bluesquares) flow solution for the x -component of the fluid phase velocity, v f , for φ = 0 . ; d) time-history of the flow solutionat x = 0 . m for several solid volume fractions; and e) comparison of analytical (dashed lines) and simulated (colored symbols)flow solutions at x = 0 . m for three sets of ∆ p f = p in − p out . To simulate this problem, we generate a regular Cartesian grid with 200 ×
10 linear elements for thenodal (FE) basis functions — with bi-linear (or “tent”; see [78]) basis functions — and discretize the poroussolid using 4,000 material points (MP) as in Figure 4b. To ensure accuracy of the fluid flow solution, wegenerate triangular finite volumes (FV) from this same Cartesian grid with additional refinement near theends of the porous solid (where ∂n/∂x is large). We then assign static pressure boundary conditions, p in and p out , at the ends of the pipe for three different pressure drops, ∆ p f ∈ { , , } kPa, and sevendifferent granular packings, φ ∈ { . , . , . , . , . , . , . } . After reaching a steady flow, as inFigure 4d, we compare the numerical value of v fx ( x ) with the analytical solution along length of the pipe.A representative flow solution for ∆ p f = 100 kPa and φ = 0 . is shown in Figure 4c, and a comparison ofthe flow solutions at x = 0 . m along the pipe for each set of pressure drops and volume fractions is shownin Figure 4e. 16 .2. Consolidation test A common numerical test for mixture simulation frameworks is the one-dimensional consolidation of afully-saturated, elastic, porous solid under external loading (see [20, 65]). Here we consider a H = 1 m tall,infinitely wide, densely packed granular column with average grain diameter d = 0 . mm and volume fraction φ = 0 . . This column is subjected to an external load, σ , of 10 kPa. The granular solid is assumed to befully saturated by water ( κ = 2 . GPa, η = 1 mPa · s, ρ f = 1000 kg/m ) and is treated as a linear elasticmaterial with Young’s modulus E ( E = 10 MPa), Poisson’s ratio ν ( ν = 0 . ), and density ρ s ( ρ s = 2650 kg/m ). The fluid is assumed to have zero pressure at z = 1 m, and the base ( z = 0 m) of the granular columnis assumed to be impermeable. Here we again use the Carman–Kozeny formula for drag.For this consolidation problem, the pore fluid pressure, p f ( z, t ) , has the following analytical form, p f ( z, t ) = ∞ (cid:88) m =0 σ M sin( M z/H ) e − M T v , with M = π m + 1) , T v = c v tH , (8.3)for c v = E v n d / (180(1 − n ) η ) and E v = E (1 − ν ) / (1 − ν − ν ) (see [65, 89]). As in the previous section,given this exact solution to the governing equations in (3.1), our numerical framework can be evaluated bycomparing its predictions to those in (8.3). Figure 5: Numerical consolidation: a) schematic diagram of test alongside finite volume (FV) grid, finite element (FE) grid,and material points (MP) used in numerical test; b) comparison of analytical (solid lines) and simulated (blue squares) porepressure solutions at T v = 0 . , . , . , . , , and . To evaluate our numerical framework as applied to this problem, we generate a 10 ×
100 element CartesianFE grid with 4 material points per grid cell (4,000 total). Standard bi-linear (or “tent”) basis functions areused in this discretization. To improve accuracy near the top surface, we generate the triangular FV meshshown in Figure 5a (with additional refinement at z = 1 m). Impermeable, smooth boundary conditions areapplied to the two walls and lower boundary of the simulated domain, and two traction boundary conditionsare imposed at the top surface: the vertical load σ is applied to the top layer of material points and a zeropressure condition is assigned to the fluid. To avoid non-physical pressure oscillations in the fluid at theonset of these loading conditions ( t = 0 ), we apply numerical damping of the form described in [20] at thevery beginning of the simulation.After simulating this problem for 2.5s of simulated time, we compare the fluid pressure solution, p f ( z, t ) ,to the analytical solution at T v = z = 1 m, the overall solution profiles are in very close agreement.17 .3. Collapse of a submerged granular column: comparison to two-phase MPM The third case that we examine using this numerical framework is the collapse of two columns of glassbeads submerged in a viscous fluid. In particular, we seek to recreate the numerical simulations reportedin [21] for two of the experiments in [90]. The two columns have the same total mass and are composedof many small, spherical glass beads ( d = 255 µ m and ρ s = 2500 kg/m ). Both columns are held behind aretaining wall and submerged in a mixture of water and Ucon oil ( η = 12 mPa · s and ρ f = 1000 kg/m ),which fills an open tank measuring 70cm × × × φ = 0 . , H = 0 . cm, and L = 6 cm and the otheris loosely packed with φ = 0 . , H = 0 . cm, and L = 6 cm. At the beginning of the experiments, theretaining walls are removed and the columns collapse under their own weight. See Figure 6a for a schematicof the initial configurations of the columns, and see Figure 6b for a schematic of their final shape along withthe location of the pressure sensor used to measure the fluid pressure, p f , at the base of the columns. -10 0 10 20 30 40-200-1000100200300 Figure 6: Collapse of a submerged column of glass beads: a) initial configuration of granular column with height H , length L , and solid volume fraction φ ; b) final configuration of granular column and position of pressure sensor; c) comparison offront position, L , found experimentally (symbols; see [90]), with MPM (dashed lines; see [21]), and with FV-MPM (solid anddash-dotted lines); d) comparison of excess fluid pressures found experimentally (symbols), with MPM (dashed lines), and withFV-MPM (solid and dash-dotted lines). * indicates FV-MPM simulation run with zero pressure BC (see Figure 7). To simulate these two experiments, we use the same elasto-plastic effective granular stress model for ˜ σ asis reported in [21], which accounts for frictional granular flow ( µ = 0 . ), rate-dependence ( µ = 1 . and b = 0 . ), shear induced dilation ( φ m = 0 . , a = 1 . , and K = 4 . ), granular elasticity ( E = 10 kPa, ν = 0 . ), and free granular separation. We model the fluid phase of the mixture as a weakly compressibleviscous fluid with bulk modulus κ = 10 kPa, baseline density ρ f = 1000 kg/m , and effective viscosity η r = η (1 + φ ) . The drag interaction is modeled using the form of ˆ F ( φ, Re) from [72].18o discretize the simulation domain, we use the 300 ×
100 element Cartesian FE grid from [21] — withthe bi-cubic-spline basis functions from [79] — and represent both granular columns with 4 material pointtracers per FE grid cell. As in [21], we model the boundaries of the domain as frictional walls. The primarydifference between our simulation approach and that used in [21] is the treatment of the fluid phase of themixture. In [21], the fluid fills the tank to a height of 8cm and is modeled using its own set of materialpoint tracers, which are shown in gray in Figure 7 and allow for direct simulation of the fluid free surface.In this work, however, we let the tank be fully filled so that the fluid phase can be modeled on a set of finitevolumes, which are defined by the same 300 ×
100 element Cartesian grid and have an impermeable upperboundary.
Time: 0s
Time: 2s
Time: 8s M P M F V - M P M F V - M P M * zero pressure surface impermeable boundaryzero pressure boundary* 𝜙 = 0.55𝜙 = 0.55 𝑎)𝜙 = 0.55 M P M F V - M P M 𝜙 = 0.60 𝜙 = 0.60 𝜙 Figure 7: Simulation snapshots for collapse of a submerged column of glass beads. The images in each column represent theflow solution at 0s, 2s, and 8s, respectively, and are found using the simulation frameworks and initial conditions labeled atthe left of the figure. The first column images indicate the relevant boundary conditions applied in each case (* indicates theapplication of a non-physical, pressure BC at the upper surface), and the arrows in the second and third columns indicatethe strength and direction of the fluid velocity field, v f . The length scale of these arrows is 0.25s × v f . The material pointsassociated with the granular phase are colored by the local granular volume fraction, φ , and the material points associated withthe fluid phase are colored gray. Inset a) highlights the regular Cartesian grid used in the MPM and FV-MPM simulations aswell as markers for the initial material point discretization. In addition to the two simulations run in this manner (labeled ‘FV-MPM’ in Figures 6 and 7), we alsoassess the potential impact of the fluid free surface on these results by running a third simulation on a30cm × ×
80 element) Cartesian grid using a zero pressure, static upper boundary (labeled ‘FV-MPM*’ in Figures 6 and 7). This zero pressure condition does not represent a physical boundary conditionbut is a closer numerical approximation to the boundary conditions used in [21]. The various boundaryconditions that are applied to the fluid phase of these mixtures are shown in the first column of Figure 7.19 comparison of the predictions made by MPM and FV-MPM with the experimental measurementsreported in [90] is shown in Figures 6c and 6d, and snapshots of these simulated flow solutions (at 0s, 2s,and 8s) are shown in Figure 7. Both simulation approaches show close agreement with the front positions(i.e. the leading edge of the collapsing columns) and excess pore pressures (i.e. the fluid pressure, p f , insidethe columns minus any hydrostatic component) reported from the experiments, and the dynamics of thecollapsing columns are remarkably similar. However, there are four differences apparent in Figures 6 and 7that should be discussed.The first two differences highlight the relative advantages of the FV-MPM approach over MPM. First,in Figure 6d the excess pore pressure found in [21] for the densely packed column is about 50 Pa lowerthan the experimental measurements and the results found in this work; this is likely an artifact of theringing instability in MPM (see [91, 92]) and is avoided by evaluating the fluid pressure with static finitevolumes. Second, the snapshots of the MPM flow solutions at 8s of simulated time in Figure 7 — after thecolumns have essentially settled — show the growth of noticeable secondary fluid flows that do not dissipateand are absent from the FV-MPM simulations; these spurious motions are indicative of quadrature errorsaccumulating in the MPM solutions (see [93]) and are examples of another numerical artifact that is avoidedwith FV-MPM.The other two differences that are apparent in these results highlight a short-coming of the proposedmethod: an inability to model the fluid free surface. First, in Figure 6c, the final front position in the‘FV-MPM’ simulations of the loosely packed column is about 2cm less than the experimental measurements,the MPM results, and the ‘FV-MPM*’ results; this is likely due the lack of a nearby zero pressure boundarycondition, which is associated with this fluid surface. Second, in Figure 7, the fluid flow solution predicted byMPM for the loosely packed column shows substantial sloshing of the mixture surface: something markedlyabsent from both the ‘FV-MPM’ and ‘FV-MPM*’ solutions. To capture the behavior of the free fluid surfacein this framework, without using Lagrangian material tracers, requires augmentation of the method withsome form of surface tracking scheme (e.g. the volume of fluid method or the level set method; see [51, 52]). We now look at problems that FV-MPM is uniquely suited to simulate and model. Consider a two-dimensional channel measuring 25cm tall and 50cm long filled with an air-like fluid ( ρ f = 1 . kg/m and ϑ f = 298 K). Suppose that a pile of sand-like material ( ρ s = 2700 kg/m , d = 0 . mm, and φ = 0 . )measuring 5cm tall and 20cm wide is situated near the front of the channel as in Figure 8a. In this numericalexample, we subject this pile of sand to a steady flow of air and model its erosion using FV-MPM.To simulate this problem, we use the identical FE and FV grids shown in Figure 8b — with linear FEbasis functions — and represent the pile of sand with 1,600 material point tracers. The sand is modeledusing the effective granular stress model for ˜ σ from [67] that accounts for simple frictional granular flow( µ = µ = 0 . ), free granular separation ( ¯ ρ c = 1620 kg/m ), and granular elasticity ( E = 16 MPa and ν = 0 . ). We treat the air-like material as an ideal gas with heat capacity ratio γ r = 1 . , specific gas constant R = 24 J/kg · K, viscosity η = 18 µ Pa · s, and coefficient of thermal conductivity k f = 0 . W/m · K. Thedrag model from [72] is again used to determine ˆ F ( φ, Re) . To capture the impact of unresolved turbulenceon our flow solution, we add the Smagorinsky eddy viscosity ( µ v = ¯ ρ f α h h α √ (cid:107) D f (cid:107) with α h ≈ . ) to ourfluid model as in [94, 95].To induce air flow through the channel, we assign pressure boundary conditions that correspond to anunimpeded air speed of 17 m/s (i.e. a stagnation pressure and temperature at the channel inlet, p ∗ f = 8700 Pa and ϑ ∗ f = 299 . K, and a static pressure and temperature at the channel outlet, p f = 8408 Pa and ϑ f = 298 K). Figure 8c and 8d show snapshots of the FV-MPM flow solution at 1.5s and 3s, respectively,including the effective air temperature, ϑ f , and equivalent granular shear strain measure, ¯ γ p , from [67]. Here ¯ γ p represents the accumulated shear deformation that the granular material has undergone. The high valueof ¯ γ p on the leeward slope of the sand pile is indicative of material that has undergone substantial shearing;in this case, it indicates that this material has been dragged along the windward face and over the crest ofthe dune. 20 ) FV/FEMP stagnation pressure inlet static pressure outletimpermeable boundary
50 cm 𝑏) 𝑐) 𝑑) 𝜗 𝑓 [°K] ҧ𝛾 𝑝 Figure 8: Numerical simulation sand pile erosion by air in two dimensions: a) schematic diagram of numerical test showingboundary conditions and initial dimensions of sand pile; b) finite volume (FV) and finite element (FE) triangular grids alongwith initial material point (MP) distribution (representing the sand pile at 0.0s); c) mixed flow solution at 1.5s with air coloredby effective temperature, ϑ f , and material points (sand) colored by the effective shear strain measure, ¯ γ p ; d) mixed flow solutionat 3.0s. Tick marks in c) and d) mark 2cm intervals. We continue the assessment of our method by extending the previous example into three dimensions.Consider a channel measuring 10cm tall, 15cm wide, and 25cm long filled with an air-like fluid ( ρ f = 1 . kg/m and ϑ f = 298 K). Suppose a conical pile of sand ( ρ s = 2700 kg/m , d = 0 . mm, and φ = 0 . )measuring 4cm tall and 12cm wide is situated near the front of the channel as in Figure 9a. In this numericalexample, we subject this pile of sand to a 1s gust of air and model its erosion using FV-MPM.To simulate this problem, we use the identical 75 × ×
30 element Cartesian FE and FV grids shown inFigure 9b — with tri-linear “tent” FE basis functions — and represent the sand pile with 32,224 materialpoint tracers. The sand is again modeled using the effective granular stress model for ˜ σ from [67] thataccounts for simple frictional granular flow ( µ = µ = 0 . ), granular separation ( ¯ ρ c = 1620 kg/m ), andgranular elasticity ( E = 16 MPa and ν = 0 . ). We treat the air-like material as an ideal gas with heatcapacity ratio γ r = 1 . , specific gas constant R = 24 J/kg · K, viscosity η = 18 µ Pa · s, and coefficient ofthermal conductivity k f = 0 . W/m · K. The drag model from [72] is used to determine ˆ F ( φ, Re) , and theSmagorinsky eddy viscosity correction from [94, 95] is applied to the fluid.The channel boundary conditions are identical to those used in the two-dimensional erosion problem andcorrespond to an unimpeded air speed of 17m/s. At the beginning of the simulation, the air is allowedto flow freely; however, after 1s of simulated time, the boundary conditions are reset to stop the flow (i.e. p ∗ f = p f = 8408 Pa and ϑ ∗ f = ϑ f = 298 K). With more degrees of freedom in three dimensions, a morecomplex fluid flow occurs. Figures 9c and 9d show snapshots of the FV-MPM flow solutions at 0.5s and 2s,respectively. In the 0.5s snapshot, streamlines are traced out beginning from the center-line of the channelinlet and colored by the local air speed, (cid:107) v f (cid:107) ; at 2s, the gust of air has stopped and there is no longerany airflow. The material points in both snapshots are colored by the accumulated equivalent shear strainmeasure, ¯ γ p . 21 ) FV/FE MP 𝑏) 𝑐) 𝑑)
25 cm stagnation pressure inlet static pressure outletimpermeable boundary 𝒇 [m/s] ҧ𝛾 𝑝 Figure 9: Numerical simulation of sand pile erosion by gust of wind in three dimensions: a) schematic diagram of numericaltest showing boundary conditions and initial dimensions of sand pile; b) finite volume (FV) and finite element (FE) Cartesiangrids along with initial material point (MP) distribution (representing the sand pile at 0.0s); c) mixed flow solution at 0.5s withair streamlines (originating along center-line of inlet) colored by velocity magnitude, || v f || , and material points (sand) coloredby the effective shear strain measure, ¯ γ p ; d) mixed flow solution at 2.0s, after wind gust has stopped. An interesting feature of this three-dimensional erosion process is the change in shape of the initiallyconical pile of sand. As the simulation progresses, the sand forms a crescent shape with two long hornsextending in the direction of the airflow. When the material has settled, the final shape of the pile isreminiscent of a barchan dune: a common dune type observed in deserts around the world. An exampleimage of several barchan dunes in Peru (from [96]) is shown in Figure 10a alongside a schematic of the classicbarchan dune features in Figure 10b. To highlight the similarity between our simulated sand pile and thischaracteristic shape, we show several of the φ = 0 . surface contours from our FV-MPM simulation in Figure10c. Although this problem is an example of erosion over a much shorter time-scale, the similarity of thesedune shapes suggests that FV-MPM is capable of modeling real world processes.22 ) 𝑏) windward side slip face horns/wings wind direction 𝑎) Figure 10: Numerical simulation of sand pile erosion by gust of wind in three dimensions: a) a field of crescent-shaped barchansand dunes in the desert between Chimbote and Casma on the coast of Peru (image source: [96]); b) schematic diagram ofcharacteristic barchan dune shape (see [97]); b) φ = 0 . surface contours for flow solutions at 0s, 0.5s, and 2s. Black contourlines in b) are marked at 4mm vertical increments, and the tick marks along the boundary denote 2.5cm increments. In this final numerical example, we examine a problem that is of growing interest to the aerospacecommunity: the effect of rocket exhaust plumes on the dusts and sediments found on extraterrestrial bodies.In the last decade many authors have investigated this effect using numerical techniques that focus on theentrained granular material (i.e. the grains that are kicked up from the granular surface; see [98, 99]). Theseapproaches generally rely on empirical erosion rates, which are based on laboratory experiments and analysesof previous extraterrestrial landings (e.g. subsonic cratering experiments [100] and analyses of the Apollo12 landing [101]), and do not account for the changing topology of the granular surface as it craters. Thissuggests that FV-MPM may be of use for modeling these types of problems, as it is capable of simulatingthe complex interaction between the impinging exhaust gases, the entrained granular material, and thedeforming, solid-like granular surface.Suppose we intend to use the Apollo Lunar Module Descent Engine (LMDE; see [102, 103]) in the designof a Martian lander and are interested in how the exhaust gases from this engine will interact with the dustyMartian surface. In particular, we want to predict the surface cratering that will occur as the lander comesto rest. In this numerical example, we simulate a rough approximation of this problem using FV-MPM.To perform this simulation, we consider the axisymmetric domain shown in the schematic in Figure 11a.The base of this domain is impermeable and covered in a 0.5m thick bed of Martian soil, and the top andoutside boundaries of this domain are open to the Martian atmosphere ( p f = 610 Pa and ϑ f = 210 K). Thesimulated LMDE is fixed in position with the nozzle exit 1.5m above the soil surface. To model the rocketexhaust accurately, we treat the nozzle wall as impermeable and friction-less, and assign a chamber pressureand temperature consistent with the operating conditions of the LMDE ( p ∗ f = 710 kPa and ϑ ∗ f = 3000 K).For our numerical implementation, we use the identical tetrahedral FE and FV wedge grids shown inFigure 11b — with linear FE basis functions — and represent the Martian soil with 498,563 material pointtracers. We model the soil using the effective granular stress model for ˜ σ from [21], which accounts forfrictional granular flow ( µ = µ = 0 . ), shear induced dilation ( φ m = 0 . , a = 1 . , and K = 4 . ),granular elasticity ( E = 10 MPa, ν = 0 . ), and free granular separation. The material properties for theMartian soil are approximated from [104], ( ρ f = 2400 kg/m and d = 0 . mm). We treat the rocket exhaust23 tmospheric conditionsimpermeable boundary axisymmetric nozzle wallcombustion conditions 𝑎) 𝑏) FV/FE MP Figure 11: Numerical simulation of rocket exhaust impinging on Martian soil: a) schematic diagram of numerical test showingboundary conditions, dimensions, and initial configuration of soil; b) finite volume (FV) and finite element (FE) tetrahedralgrids along with initial material point (MP) distribution. Simulation is performed on a 60 ◦ wedge of a cylindrical domain withsymmetric boundary conditions applied to reflected boundaries of the domain. Martian gravity (3.7 m/s ) is used for g . and Martian atmosphere as a single ideal gas with heat capacity ratio γ r = 1 . , specific gas constant R = 189 J/kg · K, viscosity η = 14 µ Pa · s, and coefficient of thermal conductivity k f = 0 . W/m · K. Thedrag model from [72] is used to determine ˆ F ( φ, Re) .There is one final consideration that must be made before we can simulate this problem in our numericalframework: it is necessary to include additional numerical dissipation to suppress oscillations near strongshocks and avoid the carbuncle phenomenon near the central axis (see [105, 106]). This requirement is notunique to our numerical method and has been studied thoroughly in the literature (e.g. see [107, 108]). Herewe extend the artificial viscosity approach from [105], which adds artificial diffusion to regions near strongshocks but does not affect smoothly varying flows. Further discussion of this approach can be found inAppendix F.With this problem correctly implemented, we begin our simulation by ramping up the chamber pressureand temperature over several milliseconds. Once the simulated engine reaches the desired operating condi-tion, we allow the flow to develop according to the equations of motion in (3.1). Figure 12 shows snapshotsof the mixed flow solution at 0.1s intervals for this rocket impingement and cratering problem. The leftcolumn of Figure 12 shows the reflections of our 60 ◦ wedge domain around the central axis and allows usto visualize the deformation of the soil bed using the material point tracers (colored by the local volumefraction, φ ). The right column of Figure 12 shows slices through the simulated domain and allows us tovisualize the shape of the exhaust plume predicted by the finite volumes (colored by the local mach number, M = (cid:107) v f (cid:107) / (cid:112) γ r Rϑ f ).Several features of this predicted flow solution are similar to those found in previous studies, includingthe mach contours inside the nozzle and the recirculation bubble just above the granular surface (see [98]);however there are several new features that have not been previously simulated. In particular, we observethat the flow direction of the exhaust gases changes by about 45 ◦ vertically as the crater is forming. Alltogether, our results highlight a strength of FV-MPM in comparison with more traditional approaches: theability to accurately simulate the coupling between surface deformations, entrained granular material, andimpinging fluids. 24 igure 12: Numerical simulation of rocket exhaust impinging on Martian soil: a) flow solution at 0.0s, b) flow solution at 0.1s,c) flow solution at 0.2s, and d) flow solution at 0.3s. Left column: renderings of Apollo Lunar Module Descent Engine (LMDE)and material point tracers colored by the local volume fraction φ . Right column: slice of simulated domain extending radiallyoutward from the center-line highlighting the overlapping fluid and granular domains; the fluid domain is colored by the localmach number M = (cid:107) v f (cid:107) / (cid:112) γ r Rϑ f . . Conclusion We have introduced a novel numerical simulation approach for a special class of engineering problems:large deformation flows of mixtures of fluids and porous solids. In this work, we have focused our attentionon fluid–sediment mixtures. These types of problems are unique in that they have features that require highaccuracy numerical advection (e.g. granular stresses, history dependent variables, porosity, etc.), suggestinga Lagrangian simulation approach, as well as features that make such Lagrangian approaches difficult (e.g.highly turbulent regions, inlet and outlet conditions, etc.).To overcome these challenges, we have proposed the FV-MPM simulation framework, which uses materialpoint tracers and an Eulerian finite element grid to solve the equations of motion of the granular phase ofthese mixtures, and an Eulerian finite volume grid to solve the equations of motion of the fluid phase of thesemixtures. By tracking the motion of the granular phase on the set of material point tracers, the advectionof masses, stresses, and material properties occurs without the numerical diffusion associated with purelyEulerian approaches. Additionally, by using finite volumes to calculate the motion of the fluid, the numericalintegration of the fluid equations is not affected by the amount of deformation that occurs over the courseof a long simulation.This numerical framework builds on recent successes using two-phase MPM approaches (see [20, 21, 29])by incorporating thoroughly studied aspects of more common FVM approaches (see [49, 55, 57]). In ouranalysis, we have shown that the proposed framework is robust and capable of addressing several classes oftime-dependent mixture problems, from soil compaction to submerged granular flows to sediment erosion byair. Additional discussion of the method, its accuracy and assumptions, and useful approximations can befound in the Appendix.
Acknowledgments
AB and KK acknowledge support from Army Research Office Grants W911NF-16-1-0440 and W911NF-19-1-0431.
AppendicesAppendix A. Mixture Theory and Simplifying Assumptions
The governing equations summarized in (3.1) derive from the mixture theories of [22, 23, 24]. In thissection, we briefly review the specialization of these theories to fluid–sediment mixtures and explicate thesimplifying assumptions that we’ve made.As stated above, this formulation assumes that the grains are quasi-mono-disperse, rough, incompressiblewith density ρ s , and fully immersed in a compressible fluid with density ρ f . A representative volumeof material, Ω , can therefore be decomposed into a solid volume, Ω s , and a fluid volume, Ω f , such that Ω = Ω s ∪ Ω f . Figure A.1 shows how this volume is decomposed and the important step of homogenizing thesolid and fluid volumes into two, overlapping continua. Appendix A.1. Homogenization
Following from the works of [21, 22, 23, 24, 75], we define the effective densities ( ¯ ρ s and ¯ ρ f ), continuumvelocities ( v s and v f ), and effective internal energies ( ε s and ε f ) such that conservation of mass, momentum,and energy in the continuum correspond to conservation of mass, momentum, and energy in the real mixture.Toward this end, we consider a representative volume of material, Ω , that contains a large number ofindividual grains: for the continuum approximation to be valid, the volume must be large enough to smooth26 igure A.1: Pictorial description of the representative volume Ω and boundary ∂ Ω , the decomposition of the domain into fluidand solid volumes, and the homogenization of the two phases. Here Ω = Ω s ∪ Ω f and ∂ Ω = ∂ Ω s ∪ ∂ Ω f with ∂ Ω ∗ defining theinterior surface separating the solid and fluid domains. out grain-scale phenomena and capture the average bulk behavior of the mixture. For such a volume, wedefine: (cid:90) Ω ¯ ρ s dv ≡ (cid:90) Ω s ρ s dv, (cid:90) Ω ¯ ρ f dv ≡ (cid:90) Ω f ρ f dv, (cid:90) Ω ¯ ρ s v s dv ≡ (cid:90) Ω s ρ s v dv, (cid:90) Ω ¯ ρ f v f dv ≡ (cid:90) Ω f ρ f v dv, (cid:90) Ω ¯ ρ s ( ε s + v s · v s ) dv ≡ (cid:90) Ω s ρ s ( ε + v · v ) dv, (cid:90) Ω ¯ ρ f ( ε f + v f · v f ) dv ≡ (cid:90) Ω f ρ f ( ε + v · v ) dv, (A.1)with ρ s and ρ f the true material densities, v the true material velocity, ε the true internal energy, and Ω s and Ω f defined as in Figure A.1. Following this procedure, the local solid volume fraction ( φ ) can be definedas the ratio of the volume of solid grains ( Ω s ) to the volume of mixture ( Ω ) and the local porosity ( n ) can bedefined similarly. If many such volumes are chosen over which variations in these local fields are negligible(i.e. relatively small volumes), we can then relate the effective fields to the true mixture densities as in (2.1). Appendix A.2. Mass Conservation
Conservation of mass within an arbitrary volume of real mixture
Ω = Ω s ∪ Ω f with boundary ∂ Ω = ∂ Ω s ∪ ∂ Ω f has the following form: ∂∂t (cid:90) Ω s ρ s dv + (cid:90) ∂ Ω s ρ s v · ˆ n da = 0 ,∂∂t (cid:90) Ω f ρ f dv + (cid:90) ∂ Ω f ρ f v · ˆ n da = 0 , (A.2)27ith ˆ n the outward facing surface normal vector, and ∂ Ω s and ∂ Ω f defining the exterior solid and fluidsurfaces as in Figure A.1. (Note that ∂ Ω ∗ , the interior surface within the domain Ω , moves with the particleboundaries, therefore there is no mass flux across it.) Using the identities in (A.1) for Ω encompassingrelatively small volumes, (A.2) can be re-expressed as, ∂∂t (cid:90) Ω ¯ ρ s dv + (cid:90) ∂ Ω ¯ ρ s v s · ˆ n da = 0 ,∂∂t (cid:90) Ω ¯ ρ f dv + (cid:90) ∂ Ω ¯ ρ f v f · ˆ n da = 0 , (A.3)and since (A.3) must hold for any representative volume, we require, ∂ ¯ ρ s ∂t + div(¯ ρ s v s ) = 0 ,∂ ¯ ρ f ∂t + div(¯ ρ f v f ) = 0 . (A.4) Appendix A.3. Momentum Conservation
Momentum conservation within an arbitrary volume of real mixture
Ω = Ω s ∪ Ω f with exterior boundary ∂ Ω = ∂ Ω s ∪ ∂ Ω f has the following form: ∂∂t (cid:90) Ω s ρ s v dv + (cid:90) ∂ Ω s ρ s v ( v · ˆ n ) da = (cid:90) Ω s ρ s g dv − (cid:90) ∂ Ω ∗ t ( ˆ n ) da + (cid:90) ∂ Ω s t ( ˆ n ) da,∂∂t (cid:90) Ω f ρ f v dv + (cid:90) ∂ Ω f ρ f v ( v · ˆ n ) da = (cid:90) Ω s ρ f g dv + (cid:90) ∂ Ω ∗ t ( ˆ n ) da + (cid:90) ∂ Ω f t ( ˆ n ) da, (A.5)with t ( ˆ n ) the surface traction vector, a function of the outward pointing surface normal vector ˆ n ; g thegravitational acceleration vector; and ∂ Ω ∗ the interior surface separating the solid and fluid domains of themixture. (Note that ∂ Ω ∗ moves with the particle boundaries, therefore there is no momentum flux acrossit.) A decomposition of the true velocity ( v ) in the solid and fluid phases of the mixture into a homogenizedcomponent ( v s and v f ) and fluctuational component ( δv s and δv f ) is common in poromechanics texts (see[22, 75]) and necessary to account for the effects of tortuosity (see [76, 77]). Here δv s ≡ v − v s in the soliddomain ( Ω s ) and δv f ≡ v − v f in the fluid domain ( Ω f ). Adding this decomposition and the definitions in(A.1) to (A.5) yields, ∂∂t (cid:90) Ω ¯ ρ s v s dv + (cid:90) ∂ Ω ¯ ρ s v s ( v s · ˆ n ) da = (cid:90) Ω ¯ ρ s g dv − (cid:90) ∂ Ω ∗ t ( ˆ n ) da + (cid:90) ∂ Ω s t ( ˆ n ) − ρ s δv s ( δv s · ˆ n ) da,∂∂t (cid:90) Ω ¯ ρ f v f dv + (cid:90) ∂ Ω ¯ ρ f v f ( v f · ˆ n ) da = (cid:90) Ω ¯ ρ f g dv + (cid:90) ∂ Ω ∗ t ( ˆ n ) da + (cid:90) ∂ Ω f t ( ˆ n ) − ρ f δv f ( δv f · ˆ n ) da. (A.6)We can evaluate the exterior surface integrals in (A.6) by defining an effective Cauchy stress for each phaseof the continuum mixture ( σ s and σ f ) according to the integral form of Cauchy’s Theorem (and including ofthe fluctuational stresses described in [22, 24]) for a sufficiently large, arbitrary domain ( Ω ) with boundary ∂ Ω = ∂ Ω s ∪ ∂ Ω f , (cid:90) ∂ Ω σ s ˆ n da ≡ (cid:90) ∂ Ω s t ( ˆ n ) − ρ s δv s ( δv s · ˆ n ) da, (cid:90) ∂ Ω σ f ˆ n da ≡ (cid:90) ∂ Ω f t ( ˆ n ) − ρ f δv f ( δv f · ˆ n ) da. (A.7)Similarly, we can evaluate the interior surface integrals in (A.6) by homogenizing the tractions along theinternal boundary ( ∂ Ω ∗ ) into a set of body forces describing the buoyant ( f b ) and drag ( f d ) interactionsoccurring between the phases, (cid:90) Ω ( f b + f d ) dv ≡ (cid:90) ∂ Ω ∗ t ( ˆ n ) da. (A.8)28ogether (A.6), (A.7), and (A.8) must hold for any volume Ω such that momentum conservation can beexpressed locally as, ∂ ¯ ρ s v s ∂t + div (cid:0) ¯ ρ s v s ⊗ v s (cid:1) = ¯ ρ s g − f b − f d + div σ s ,∂ ¯ ρ f v f ∂t + div (cid:0) ¯ ρ f v f ⊗ v f (cid:1) = ¯ ρ f g + f b + f d + div σ f . (A.9) Appendix A.4. Energy Conservation
Conservation of energy over an arbitrary volume of real mixture
Ω = Ω s ∪ Ω f with boundary ∂ Ω = ∂ Ω s ∪ ∂ Ω f has the following form: ∂∂t (cid:90) Ω s ρ s (cid:0) ε + v · v (cid:1) dv = − (cid:90) ∂ Ω s ρ s (cid:0) ε + v · v (cid:1) ( v · ˆ n ) da − (cid:90) ∂ Ω s ( q · ˆ n ) da + (cid:90) ∂ Ω ∗ ( q · ˆ n ) da + (cid:90) Ω s q dv + (cid:90) ∂ Ω s ( t ( ˆ n ) · v ) da − (cid:90) ∂ Ω ∗ ( t ( ˆ n ) · v ) da + (cid:90) Ω s ( ρ s g · v ) dv, (A.10)and, ∂∂t (cid:90) Ω f ρ f (cid:0) ε + v · v (cid:1) dv = − (cid:90) ∂ Ω f ρ f (cid:0) ε + v · v (cid:1) ( v · ˆ n ) da − (cid:90) ∂ Ω f ( q · ˆ n ) da − (cid:90) ∂ Ω ∗ ( q · ˆ n ) da + (cid:90) Ω f q dv + (cid:90) ∂ Ω f ( t ( ˆ n ) · v ) da + (cid:90) ∂ Ω ∗ ( t ( ˆ n ) · v ) da + (cid:90) Ω f ( ρ f g · v ) dv, (A.11)where q is the true heat flow vector and q is the true rate of internal heat generation. (Note that althoughthere is no energy flux across ∂ Ω ∗ , heat flow through this boundary is still possible.) To evaluate the exteriorboundary integrals in (A.10) and (A.11), we first define a set of effective heat flow vectors ( q s and q f ) thatinclude the extra mixing of internal energies due to local fluctuations in velocity as, (cid:90) ∂ Ω ( q s · ˆ n ) da ≡ (cid:90) ∂ Ω s ( q · ˆ n ) + ρ s ε ( δv s · ˆ n ) da, (cid:90) ∂ Ω ( q f · ˆ n ) da ≡ (cid:90) ∂ Ω f ( q · ˆ n ) + ρ f ε ( δv f · ˆ n ) da, (A.12)and a set of effective phase-wise rates of internal heat generation ( q s and q f ) as, (cid:90) Ω q s dv ≡ (cid:90) Ω s q dv, (cid:90) Ω q f dv ≡ (cid:90) Ω f q dv. (A.13)29imilarly, we can evaluate the internal boundary integrals in (A.10) and (A.11) by defining an effectiveinter-phase heat flow rate ( q i ) that represents the flow from the solid domain to the fluid domain as, (cid:90) Ω q i dv ≡ − (cid:90) ∂ Ω ∗ ( q · ˆ n ) da. (A.14)Substituting these definitions along with those from (A.1), (A.7), and (A.8) into (A.10) and (A.11) allowsus to express the evolution of the effective phase-wise energies as follows, ∂∂t (cid:90) Ω ¯ ρ s (cid:0) ε s + v s · v s (cid:1) dv = − (cid:90) ∂ Ω ¯ ρ s (cid:0) ε s + v s · v s (cid:1) ( v s · ˆ n ) da + (cid:90) ∂ Ω ( σ s v s − q s ) · ˆ n da + (cid:90) Ω q s − q i + (¯ ρ s g ) · v s − ( f b + f d ) · v s dv + (cid:90) ∂ Ω s ( t ( ˆ n ) · δv s ) + ρ s ( δv s · δv s )( δv s · ˆ n ) da − (cid:90) ∂ Ω ∗ ( t ( ˆ n ) · δv s ) da, (A.15)and, ∂∂t (cid:90) Ω ¯ ρ f (cid:0) ε f + v f · v f (cid:1) dv = − (cid:90) ∂ Ω ¯ ρ f (cid:0) ε f + v f · v f (cid:1) ( v f · ˆ n ) da + (cid:90) ∂ Ω ( σ f v f − q f ) · ˆ n da + (cid:90) Ω q f + q i + (¯ ρ f g ) · v f + ( f b + f d ) · v s dv + (cid:90) ∂ Ω f ( t ( ˆ n ) · δv f ) + ρ f ( δv f · δv f )( δv f · ˆ n ) da + (cid:90) ∂ Ω ∗ ( t ( ˆ n ) · δv s ) da. (A.16)The remaining integrals over the boundaries ∂ Ω s , ∂ Ω f , and ∂ Ω ∗ in (A.15) and (A.16) capture the extradissipation and diffusion of energy that occurs to do the fluctuations in the true velocity field. To accountfor this in our equations of motion, we introduce measures of this extra work in the solid phase ( δw s ), in thefluid phase ( δw f ), and between the phases ( δw i ) as follows, (cid:90) Ω δw s ≡ (cid:90) ∂ Ω s ( t ( ˆ n ) · δv s ) + ρ s ( δv s · δv s )( δv s · ˆ n ) da, (cid:90) Ω δw f ≡ (cid:90) ∂ Ω f ( t ( ˆ n ) · δv f ) + ρ f ( δv f · δv f )( δv f · ˆ n ) da, (cid:90) Ω δw i ≡ (cid:90) ∂ Ω ∗ ( t ( ˆ n ) · δv s ) da. (A.17)For large domains Ω over which δv s and δv f have no inherent directionality, these integrals will vanish;however, if the local fluctuations in velocity have directionality that is correlated with their magnitudes orwith the local traction vector, t ( ˆ n ) , these extra contributions to the evolution of internal energy must betaken into account. Together (A.15), (A.16), and (A.17) must hold for any volume Ω such that energyconservation can be expressed locally as, ∂ ¯ ρ s E s ∂t = − div(¯ ρ s E s v s ) + div( σ s v s − q s ) + q s − q i + (¯ ρ s g ) · v s − ( f b + f d ) · v s + δw s − δw i ,∂ ¯ ρ f E f ∂t = − div(¯ ρ f E f v f ) + div( σ f v f − q f ) + q f + q i + (¯ ρ f g ) · v f + ( f b + f d ) · v s + δw f + δw i , (A.18)30ith E s ≡ ε s + v s · v s and E f ≡ ε f + v f · v f the specific total energies for each phase. Appendix A.5. Common Expressions for the Cauchy Stresses and Buoyant Force
The effective Cauchy stresses ( σ s and σ f ) defined in (A.7) are evaluated over many large domains Ω (or,as in [75], averages over an ensemble of realizations), we therefore expect these stresses contain contributionsfrom both phases of the real mixture and can be decomposed as follows. The solid phase stress, σ s , takesthe following form: σ s = ˜ σ − (1 − n ) p f I , (A.19)with ˜ σ the portion of the solid phase stress resulting from true granular contacts and microscopic viscousinteractions between grains due to immersion in a fluid medium (e.g. lubrication forces), and p f the true fluid pore pressure. The fluid phase stress, σ f , takes the following form: σ f = τ f − np f I , (A.20)with τ f a deviatoric shear stress tensor and np f I the spherical stress given by the porosity of the mixtureand true fluid pore pressure.Additionally, we decompose the interior traction integral in (A.8) into normal components (i.e. buoyantforces; f b ) and shear components (i.e. drag forces; f d ). Following the works of [23, 24], we let the buoyantforce take the following form: f b = p f ∇ n, (A.21)which accounts for the distribution of pressures on particle surfaces that arises due to the presence of a porefluid. The drag force ( f d ) is often defined empirically as in (2.6). Appendix A.6. Closure of Equations and Simplifying Assumptions
As stated in [24], closure of this system of conservation equations requires definitions of following quan-tities in terms of the homogenized fields in (A.1): ˜ σ , p f , τ f , f d , q s , q f , q s , q f , q i , δw s , δw f , δw i . (A.22)In the main text, we have reduced this set by making the following simplifying assumptions.Since we are primarily interested in fluid–sediment flows where the absolute temperature of the grainsdoes not considerably affect the solution, we model the granular phase as an isothermal solid: the granulartemperature, ϑ s , is constant; the solid phase heat generation and heat flux, q s and q s , are implicitly defined;and the effective granular stress ˜ σ evolves according to an admissible isothermal equation of state (see [21]).If we further assume that there is no inter-phase heat flow, q i → , (A.23)then there is no need to explicitly track ε s in the mixture.In addition to modeling the granular phase as isothermal, we will also limit ourselves to consideringchemically stable mixtures ( q f → ) and neglect the extra work terms ( δw s , δw f , and δw i ), which derivefrom higher order correlations of velocity fluctuations δv s and δv f . Although there are constitutive modelsfor the average magnitudes of δv s and δv f (see [77, 109, 110]), there are few descriptions of how the directionand magnitude of these local velocity fluctuations are correlated.With these assumptions, (A.4), (A.9), and (A.18) can be simplified to the system of equations presentedin (3.1). Appendix B. Useful Numerical Approximations
In this section, we present numerical approximations that we have found useful when implementing thenumerical algorithm described in the main text. These approximations of the buoyant force vectors ( f buoy i and F buoy α ), the drag force vectors ( f drag i and F drag α ), and the fluid phase gradients ( (cid:104)∇ ρ f (cid:105) α , (cid:104)∇ ρ f v f (cid:105) α ,and (cid:104)∇ ρ f E f (cid:105) α ) improve the stability of the numerical algorithm without adding significant approximationerrors. 31 ppendix B.1. Buoyant Force Vectors The buoyant force vectors, f buoy i and F buoy α , are presented in (3.16) and (3.17), respectively. Althoughthe expression for F buoy α can be evaluated explicitly using numerical quadrature (i.e. sampling the solid andfluid phase fields within each grid cell), the exact expression for f buoy i includes a gradient that can be trickyto evaluate on non-uniform grids: f buoy i = N n (cid:88) j =1 (1 − n j ) (cid:90) Ω ∇ (cid:0) N i ( x ) N j ( x ) (cid:1) p f dv. Additionally, since this expression reflects an exact evaluation of the weak-form equation in (3.3), it mayproduce instabilities in the method when the diagonal matrices, [ M D ] and [ B D ] , are used in place of theconsistent matrices, [ M ] and [ B ] , (see [111]).To avoid these issues, we introduce the following approximation of f buoy i : f buoy i ∗ = (1 − n i ) (cid:90) Ω ∇N i ( x ) p f dv. (B.1)Here, the integral expression is easily evaluated with numerical quadrature, and we can ensure that thecoupling force is numerically stable (i.e. f buoy i ∗ = when n i = 1 ). If we assume that for a particular choiceof finite element basis functions {N i ( x ) } there exists a unique set of coefficients { ψ i } that produce the nodevalues { ψ ( x i ) } (with x i the position of the i th node) and that these functions are only non-zero over a limitedregion of the simulated domain Ω (e.g. “tent” functions or cubic splines; see [78, 79]), then it is possible toshow that, (cid:107) f buoy i − f buoy i ∗ (cid:107) ≤ c h i (cid:107)∇ n (cid:107) ∞ (cid:107)∇ p f (cid:107) ∞ (cid:90) Ω N i ( x ) dv, (B.2)with c a basis specific constant, h i the characteristic length of the finite element grid (e.g. the averagedistance between grid nodes), and (cid:107) ψ (cid:107) ∞ the maximum component of ψ across the entire simulated domain. Appendix B.2. Drag Force Vectors
The drag force vectors, f drag i and F drag α , are presented in (3.16) and (3.17), respectively. Although bothexpressions can be evaluated using numerical quadrature, we find that this approach can develop suddeninstabilities. One source of these instabilities is the consistency of the exact expression for f drag i with theweak-form equation in (3.3): in numerical algorithms that use the diagonal matrices, [ M D ] and [ B D ] , exactcalculation of these force vectors can produce spurious accelerations.To help avoid these instabilities we approximate these drag force vectors by defining a resistance field, K , as follows: K ≡ (cid:107) f d (cid:107)(cid:107) v s − v f (cid:107) . (B.3)(Note that this definition of resistance is only valid when the drag model for f d is ‘co-directional’ with thevelocity difference, v s − v f .) We then project this field onto the finite element basis and call this projection K ∗ : K ∗ = N n (cid:88) i =1 K ∗ i ( t ) N i ( x ) . (B.4)We let the coefficient K ∗ i be approximated as, K ∗ i ≈ n i (1 − n i ) η d ˆ F (1 − n i , Re ∗ i ) , (B.5)32ith Re ∗ i an approximation of the Reynolds number ( Re ∗ i ≡ (¯ ρ ∗ fi d (cid:107) v si − v ∗ fi (cid:107) ) /η ), ¯ ρ ∗ fi an approximationof the effective fluid density ( ¯ ρ ∗ fi ≡ (cid:80) N v α =1 ( V α A iα (cid:104) ¯ ρ f (cid:105) α ) /V i ), v ∗ fi an approximation of the fluid velocity( v ∗ fi ≡ (cid:80) N v α =1 ( V α A iα (cid:104) ¯ ρ f v f (cid:105) α / (cid:104) ¯ ρ f (cid:105) α ) /V i ), and V i = (cid:90) Ω N i ( x ) dv. The drag vectors can then be approximated as, f drag i ∗ = − V i K ∗ i ( v si − v ∗ fi ) , (B.6)and, F drag α ∗ = N n (cid:88) i =1 A iα K ∗ i (cid:0) v si − (cid:104) ¯ ρ f v f (cid:105) α / (cid:104) ¯ ρ f (cid:105) α (cid:1) · v si . (B.7)If we assume that for a particular choice of finite element basis functions {N i ( x ) } there exists a uniqueset of coefficients { ψ i } that produce the node values { ψ ( x i ) } and that these functions are only non-zero overa limited region of the simulated domain Ω , then it is possible to show that, (cid:107) f drag i − f drag i ∗ (cid:107) ≤ c h i (cid:107)∇ f d (cid:107) ∞ V i + c h α (cid:107)∇ v f (cid:107) ∞ (cid:107) ∂ f d /∂ v f (cid:107) ∞ V i + c h α (cid:107)∇ ¯ ρ f (cid:107) ∞ (cid:107) ∂ f d /∂ ¯ ρ f (cid:107) ∞ V i , (B.8)and, (cid:107) F drag α − F drag α ∗ (cid:107) ≤ c h i (cid:107)∇ ( K v s ) (cid:107) ∞ + c h α (cid:107) ∂ v f /∂ x (cid:107) ∞ (cid:107) ∂ f d /∂ v f (cid:107) ∞ + c h α (cid:107) ∂ ¯ ρ f /∂ x (cid:107) ∞ (cid:107) ∂ f d /∂ ¯ ρ f (cid:107) ∞ + c h i h α (cid:107)∇ v f (cid:107) ∞ (cid:107) ∂ f d /∂ v f (cid:107) ∞ + c h i h α (cid:107)∇ ¯ ρ f (cid:107) ∞ (cid:107) ∂ f d /∂ ¯ ρ f (cid:107) ∞ , (B.9)with c , c , c , and c , grid specific constants; h i the characteristic length of the finite element grid; and h α the characteristic length of the finite volume grid. If a finite element basis is chosen that has the property N i ( x i ) = 1 , then the leading error term above becomes c l h li (cid:107) ∂ ( l +1) ( K v s ) /∂ x ( l +1) (cid:107) ∞ with l the order of thebasis functions. Appendix B.3. Fluid Phase Gradients
The numerical reconstruction of the fluid phase fields ( ¯ ρ f , ρ f , ρ f v f , and ρ f E f ) shown in (4.3) uses localapproximations of the fluid phase gradients ( (cid:104)∇ ρ f (cid:105) α , (cid:104)∇ ρ f v f (cid:105) α , and (cid:104)∇ ρ f E f (cid:105) α ) computed from surroundingcentroid data. Although these reconstructions are accurate to second order, such an approach can producesignificant errors when the porosity field has unresolved fluctuations (i.e. fluctuations with a characteristiclength that is the same magnitude as the finite volume length scale, h α ). In this section, we propose analternative approach for gradient calculations which will have higher accuracy in steady, adiabatic flows,but may introduce small errors in transient problems. This alternative approach is used in the simulationsreported in sections 8.1, 8.2, 8.4, and 8.5.Consider, the one-dimensional finite volume grid and porosity field shown in Figure B.1. In steady flow,mass conservation dictates that the true momentum field must fluctuate with the porosity field, which herecorresponds to linear flow area (see the solid red line in the figure); as the porosity decreases, there shouldbe an inversely proportional change in true fluid velocity. This relationship can be expressed for adiabatic,33nviscid flows as follows: ∂ρ f ∂s = ρ f n (cid:18) M − M (cid:19) ∂n∂s ,∂ρ f v f ∂s = − ρ f v f n ∂n∂s ,∂ρ f E f ∂s = ρ f E f + p f − ρ f a n (cid:18) M − M (cid:19) ∂n∂s , (B.10)with M the local Mach number, a the local speed of sound in the fluid, and ∂s denoting infinitesimalincrements along a streamline. Figure B.1: Pictorial representation of a one-dimensional, steady-state flow solution. The blue solid line shows a representativeporosity field, n , with fluctuations that are not resolved on the finite volume grid, and the blue ‘ × ’s denote the cell-wise averagevalues. The red solid line shows a steady-state true momentum field ( ρ f v f ) solution that is consistent with the porosity fieldbelow, and the red ‘ (cid:3) ’s denote the cell-wise average values. The red long-dashed line shows the reconstructed true momentumfield found using the gradient approximations from (4.4). The red short-dashed line shows the reconstructed true momentumfield found using the corrected gradient approximations from (B.15). Applying the reconstruction procedure described in the main text to this one-dimensional problem wouldproduce a uniform momentum field, inconsistent with the grid-scale fluctuations expected in the true solution(see the long-dashed red line in Figure B.1). Although refining the finite volume grid would improve theaccuracy of the reconstruction, it is also possible to directly incorporate the analytical expressions from(B.10) into our reconstruction procedure. To do this, we calculate the cell-wise average porosity gradient, (cid:104)∇ n (cid:105) α , as follows, (cid:104)∇ n (cid:105) α = 1 V α N i (cid:88) i =1 n i (cid:90) Ω α ∇N i ( x ) dv, (B.11)and calculate the associated best estimate gradient [ ∇ n ] α and flux limiter Φ n from the surrounding centroiddata. We can then quantify the amount of grid-scale porosity fluctuation with δ ∇ n , the porosity gradientcorrection vector: δ ∇ n = (cid:104)∇ n (cid:105) α − Φ n [ ∇ n ] α . (B.12)Next, we estimate the relative motion of the solid and fluid phases by calculating the average solid phase34elocity in the finite volume cell, (cid:104) v s (cid:105) α , as, (cid:104) v s (cid:105) α = N i (cid:88) i =1 A iα v si , (B.13)and defining v ∗ , the relative fluid velocity vector, as follows, v ∗ = (cid:104) ρ f v f (cid:105) α (cid:104) ρ f (cid:105) α − (cid:104) v s (cid:105) α . (B.14)For subsonic flows, we then apply the following correction to (4.4): (cid:104)∇ ρ f (cid:105) α = Φ ρ f [ ∇ ρ f ] α + (cid:104) ρ f (cid:105) α (cid:104) n (cid:105) α (cid:18) M ∗ − M ∗ (cid:19) v ∗ ( v ∗ · δ ∇ n ) (cid:107) v ∗ (cid:107) , (cid:104)∇ ρ f v f (cid:105) α = Φ v f [ ∇ ρ f v f ] α − (cid:104) ρ f v f (cid:105) α (cid:104) n (cid:105) α (cid:18) v ∗ ( v ∗ · δ ∇ n ) (cid:107) v ∗ (cid:107) (cid:19) , (cid:104)∇ ρ f E f (cid:105) α = Φ E f [ ∇ ρ f E f ] α + ( (cid:104) ρ f E f (cid:105) α + p f − (cid:104) ρ f (cid:105) α a ) (cid:104) n (cid:105) α (cid:18) M ∗ − M ∗ (cid:19) v ∗ ( v ∗ · δ ∇ n ) (cid:107) v ∗ (cid:107) , (B.15)with p f and a calculated at the centroid of Ω α and M ∗ the relative fluid mach number, M ∗ = (cid:107) v ∗ (cid:107) /a .Applying the corrected reconstruction from (B.15) for the steady flow problem described in Figure B.1produces the field shown by the short-dashed red line in the figure: a much closer match to the analyticalsolution. For some fluid flows, particularly those near steady-state, this corrected reconstruction approachcan improve the overall accuracy of the method without requiring additional refinement. Further, for generalflows, any errors introduced by applying this correction will disappear when the porosity field is properlyresolved: (cid:13)(cid:13) ∇ ρ f (cid:12)(cid:12) X α − (cid:104)∇ ρ f (cid:105) α (cid:13)(cid:13) ≤ c h α (cid:107) ∂ ρ f /∂ x (cid:107) ∞ + A (cid:0) c h α (cid:107) ∂ n/∂ x (cid:107) ∞ + c h α (cid:107) ∂ n/∂ x (cid:107) ∞ (cid:1) , (cid:13)(cid:13) ∇ ρ f v f (cid:12)(cid:12) X α − (cid:104)∇ ρ f v f (cid:105) α (cid:13)(cid:13) ≤ c h α (cid:107) ∂ ( ρ f v f ) /∂ x (cid:107) ∞ + A (cid:0) c h α (cid:107) ∂ n/∂ x (cid:107) ∞ + c h α (cid:107) ∂ n/∂ x (cid:107) ∞ (cid:1) , (cid:13)(cid:13) ∇ ρ f E f (cid:12)(cid:12) X α − (cid:104)∇ ρ f E f (cid:105) α (cid:13)(cid:13) ≤ c h α (cid:107) ∂ ( ρ f E f ) /∂ x (cid:107) ∞ + A (cid:0) c h α (cid:107) ∂ n/∂ x (cid:107) ∞ + c h α (cid:107) ∂ n/∂ x (cid:107) ∞ (cid:1) . (B.16)Here c and c represent grid specific constants, h α represents the characteristic finite volume length scale,and { A j } represent the leading coefficients in (B.15). Appendix C. Notes on Numerical Consistency
In this section, we briefly discuss the accuracy of the numerical method described in the main text. Webegin by evaluating the consistency of solutions to the weak- and integral-forms of the governing equations in(3.2), (3.3), and (3.4) with solutions to the strong forms in (3.1). This analysis provides some basic intuitionabout the expected convergence rates of the numerical algorithm, which we then validate using the methodof manufactured solutions (MMS; see [39, 112]).
Appendix C.1. Order of Accuracy
We begin this analysis by considering the weak expression for solid phase mass conservation in (3.2).Using the definition of w from (3.8) and introducing the exact solution fields ¯ ρ † s , v † s , n † , which solve (3.1),we can express the rate of growth of density error as follows, (cid:90) Ω p D s (¯ ρ † s − ¯ ρ s ) Dt dv = − (cid:90) Ω p (¯ ρ † s − ¯ ρ s ) div( v † s ) + ¯ ρ s div( v † s − v s ) + ( v † s − v s ) · ∇ ¯ ρ † s dv, (C.1)35ith Ω p the domain of the p th material point characteristic function. If we let e ¯ ρ s ≡ ¯ ρ † s − ¯ ρ s and e v s ≡ v † s − v s it is possible to show that this growth rate is bounded as, ddt (cid:107) e ¯ ρ s (cid:107) ∞ ≤ (cid:107) e ¯ ρ s (cid:107) ∞ (cid:107)∇ v † s (cid:107) ∞ + (cid:107) ¯ ρ s (cid:107) ∞ (cid:0) c h − i (cid:107) e v s (cid:107) ∞ + c h li (cid:107) ∂ ( l +1) v † s /∂ x ( l +1) (cid:107) ∞ (cid:1) + ρ s (cid:107) e v s (cid:107) ∞ (cid:107)∇ n † (cid:107) ∞ + c h p ρ s (cid:0) (cid:107)∇ n † (cid:107) ∞ (cid:107)∇ v s (cid:107) ∞ + (cid:107)∇ ∂n † /∂t (cid:107) ∞ (cid:1) , (C.2)with c , c , and c discretization specific constants, h p the characteristic length scale of material points, h i the characteristic length of the finite element grid, and l the order of the finite element basis functions. Ifwe assume that e v s only accounts for discretization error (i.e. the error between the exact solution and theprojection of the exact solution onto the finite element basis), and that e ¯ ρ s is initially zero, then we find thatas the finite element grid and material points are refined, e ¯ ρ s will develop errors that are O ( h li ) + O ( h p ) .Now consider the weak expression for solid phase momentum conservation in (3.3). Using the definitionof w from (3.10) and introducing the exact solution fields ˜ σ † , f d † , and p † f , which solve (3.1), we can expressthe rate of growth of velocity error as follows: (cid:90) Ω ¯ ρ † s w · D s ( v † s − v s ) Dt dv = (cid:90) Ω − (¯ ρ † s − ¯ ρ s ) a s · w − ( v † s − v s ) · ( ∇ v † s ) w − ( ˜ σ † − ˜ σ ) : ∇ w + (¯ ρ † s − ¯ ρ s ) g · w − ( f d † − f d ) · w s − ( p † f − p f ) div( n † w ) − p f div (cid:0) ( n † − n ) w (cid:1) dv. (C.3)Following a similar procedure to that described above, it is possible to convert (C.3) into an inequality thatrelates the growth of solid velocity error ( e v s ) to the finite element grid scale ( h i ); the material point lengthscale ( h p ); the maximum material stretch λ max ; the solution gradients (e.g. ∇ n † and ∇ v † s ); and the errorspresent in the density, stress, drag, and pressure calculations. Under ideal initial conditions, e v s can beshown to develop errors that are O ( h li ) + O ( h p /h i ) + O ( h α /h i ) , with h α the characteristic length of the finitevolume grid.Finally, we consider the integral expressions for fluid phase mass, momentum, and energy conservationin (3.4). Introducing the exact solution fields ¯ ρ † f , v † f , E † f , τ f † , and q † f , which solve the strong form of thegoverning equations in (3.1), we can express the rates of change of the local density error as, ∂∂t (cid:18) (cid:90) Ω α ¯ ρ † f − ¯ ρ f dv (cid:19) = − (cid:90) ∂ Ω α (¯ ρ † f v † f − ¯ ρ f v f ) · ˆ n da, (C.4)the local momentum error as, ∂∂t (cid:18) (cid:90) Ω α (¯ ρ † f v † f − ¯ ρ f v f ) dv (cid:19) = − (cid:90) ∂ Ω α (cid:0) (¯ ρ † f v † f − ¯ ρ f v f )( v f · ˆ n ) + ¯ ρ † f v † f ( v † f − v f ) · ˆ n (cid:1) da + (cid:90) ∂ Ω α ( τ f † − τ f ) ˆ n − ( n † − n ) p f ˆ n − n † ( p † f − p f ) ˆ n da + (cid:90) Ω α (cid:0) (¯ ρ † f − ¯ ρ f ) g + ( f d † − f d ) + ( p † f − p f ) ∇ n + p † f ( ∇ n † − ∇ n ) (cid:1) dv, (C.5)36nd the local energy error as, ∂∂t (cid:18) (cid:90) Ω α (¯ ρ † f E † f − ¯ ρ f E f ) dv (cid:19) = − (cid:90) ∂ Ω α (cid:0) (¯ ρ † f E † f − ¯ ρ f E f ) + ( n † − n ) p f + n † ( p † f − p f ) (cid:1) ( v f · ˆ n ) da − (cid:90) ∂ Ω α (¯ ρ † f E † f + n † p † f )( v † f − v f ) · ˆ n da + (cid:90) ∂ Ω α (cid:0) ( τ f † − τ f ) v f + τ f † ( v † f − v f ) − ( q † f − q f ) (cid:1) · ˆ n da + (cid:90) Ω α (¯ ρ † f v † f − ¯ ρ f v f ) · g dv + (cid:90) Ω α (cid:0) ( p † f − p f ) ∇ n + p † f ( ∇ n † − ∇ n ) + ( f d † − f d ) (cid:1) · v † s dv + (cid:90) Ω α (cid:0) ( p † f − p f ) ∇ n + p † f ( ∇ n † − ∇ n ) + ( f d † − f d ) (cid:1) · ( v † s − v s ) dv. (C.6)Again following the procedure used to calculate C.2, it is possible to convert (C.4), (C.5), and (C.6) into aninequality that relates the growth of the fluid density error ( e ¯ ρ f ≡ ¯ ρ † f − ¯ ρ f ), velocity error ( e v f ≡ v † f − v f ),and energy error ( e E f ≡ E † f − E f ) to the finite element grid scale ( h i ); the solution gradients (e.g. ∇ n † ); andthe errors present in the density, stress, drag, pressure, and heat flux calculations. Under ideal conditions,the fluid fields can be shown to develop errors that are O ( h li ) + O ( h α ) .All together, this simplified analysis provides insight into the potential sources of error and the conver-gence rates we might expect when the method is implemented. It is clear that the initial discretization errorsin ¯ ρ s , v s , ¯ ρ f , v f , and E f — as well as the propagation of these errors to the derived quantities n , ˜ σ , p f , τ f , f d , and q f — play a crucial role determining the overall accuracy of the method. Additionally, error in any of these terms will contribute to the growth of the error of all of these terms. Therefore, assuming that anideal initial discretization is implemented and that the weak- and integral-forms of the governing equationsin (3.2), (3.3), and (3.4) are evaluated exactly, the analysis above indicates that the leading error terms inour numerical solutions should be of the following order: O ( h li ) + O ( h α ) + O ( h p ) + O ( h p /h i ) + O ( h α /h i ) .As mentioned in the main text, some of the numerical approximations that we use will introduce additionalerrors and reduce the order of accuracy of our method. Several of these approximations have been studiedthoroughly in the literature (e.g. uGIMP integration; see [39, 78]) and a few were discussed in the previoussection. With these approximations incorporated, the overall accuracy of the method should be of thefollowing order: O ( h i ) + O ( h α ) + O ( h p ) + O ( h p /h i ) . Appendix C.2. Validation with Method of Manufactured Solutions
In this section, we validate the convergence rates described above using the method of manufacturedsolutions (MMS). In the MMS, a solution to the governing equations is assumed a priori ; then the set ofexternal forces and heats necessary to achieve this solution are determined analytically (see [39, 112]). Forthe purposes of this numerical test, we adjust the governing equations in (3.1) as follows: D s ¯ ρ s Dt = − ¯ ρ s div( v s ) , ¯ ρ s D s v s Dt = div( ˜ σ ) + ¯ ρ s b s − f d − φ ∇ p f ,∂ ¯ ρ f ∂t = − div(¯ ρ f v f ) ,∂ ¯ ρ f v f ∂t = − div (cid:0) ¯ ρ f v f ⊗ v f + np f I (cid:1) + div( τ f ) + ¯ ρ f b f + f d + p f ∇ n,∂ ¯ ρ f E f ∂t = − div (cid:0) (¯ ρ f E f + np f ) v f (cid:1) + div( τ f v f ) + (¯ ρ f b f ) · v f + (cid:0) p f ∇ n + f d (cid:1) · v s + q f . (C.7)37ere b s and b f represent the solid and fluid phase external force vectors, respectively, and q f representsthe external heat flow into the fluid phase. We let the effective granular stress obey a pseudo-Neo-Hookeanmaterial model, ˜ σ := G B + K log( J ) I , (C.8)with B the deviator of the left Cauchy–Green tensor ( B = F F (cid:62) for the deformation gradient F ), J ameasure of volume change ( J = det( B ) ), and G and K the shear and bulk moduli ( G = 3 . kPa and K = 8 . kPa); additionally, we let the fluid phase obey a standard ideal gas law, p f := ρ f Rϑ f , ρ f ε f := p f γ r − , and τ f := η (1 + φ ) (cid:0) ∇ v f + ∇ v (cid:62) f − div( v f ) I (cid:1) , (C.9)with R the specific gas constant ( R = 0 . J/kg · K), γ r the ratio of specific heats ( γ r = 1 . ), and eta the dynamic fluid viscosity ( η = 2 Pa · s). The final component of our mixture model is the inter-phase drag, f d , which we let take the form described in [72] for d = 0 . m.To calculate the MMS, we first choose the desired solution fields for the solid and fluid phases of themixture. For simplicity, we choose uniform density fields, ρ † s = 1000 kg/m and ρ † f = 117 . kg/m ; auniform porosity field, n † = 0 . ; and a uniform internal fluid energy field, ε † f = 21 . J/kg. We then choosetwo independent, divergence-free velocity fields: v † f = A ( t ) ( x − y ( y − − x ( x − y − , A ( t ) = − t , (C.10)and, v † s = B ( t ) − y ( x + y − x ( x + y − if x + y ≤ , B ( t ) = 4 t , (C.11)with v † s = for x + y > . These velocity fields also have the useful property that v † s = v † f = for x = ± or y = ± . Calculation of the necessary external forces ( b s and b f ) and heat flow ( q f ) to produce thesesolutions can be found by carefully inverting the system of equations in (C.7) and calculating the deformationgradient as follows, F = ∂ x ∂ X , and D s x Dt = v s , (C.12)with x a spatial point in the deformed solid phase material, which corresponds to a position X in thereference (or initial) material configuration.With the analytical external forces and heat determined, we implement this problem in our FV-MPMframework. We use numerical quadrature to integrate the fluid body force and heat addition and — asis standard in MPM — apply the solid phase body force to the material points directly. This problemis simulated on a 2m ×
2m domain centered at the origin and discretized using three overlapping, regular,Cartesian grids. The respective grid spacing is h p for the material points, h i for the finite element grid,and h α for the finite volume grid. Standard bi-linear (or “tent”) functions are used to define the finiteelement basis, N i ( x ) . Additionally, to enable larger time-increments ( ∆ t ), we use the alternative numericalalgorithm detailed in Appendix E to integrate the time derivatives of each field. The simulated domain,boundary conditions, and results are shown in Figure C.1.To analyze the convergence behavior of our method, we simulate this problem on twelve unique sets ofmeshes and calculate the L norm of the velocity errors ( e v s and e v f ) for each as follows, (cid:107) e v s (cid:107) L = (cid:20) (cid:90) Ω (cid:107) v † s − v s (cid:107) dv (cid:21) , and (cid:107) e v f (cid:107) L = (cid:20) (cid:90) Ω (cid:107) v † f − v f (cid:107) dv (cid:21) . (C.13)The twelve sets of meshes are divided into three groups, one to test each of the convergence rates expectedof our method: O ( h i ) + O ( h α ) + O ( h p ) + O ( h p /h i ) . (Note that the latter two convergence rates are grouped38 𝑠 = 𝒗 𝑓 = 𝟎 𝑎) 𝑏) FV FE MP ℎ 𝛼 ℎ 𝑖 ℎ 𝑝 𝑐) 𝑑) 𝑒) Figure C.1: Method of manufactured solutions (MMS) results. a) Simulated domain and boundary conditions. b) FV, FE,and MP discretizations and associated characteristic lengths h α , h i , and h p . c) L -error vs. finite volume length scale, h α ; h i and h p are constant. d) L -error vs. finite element length scale, h i ; h α and h p /h i are constant. e) L -error vs. material pointrefinement measure, h p /h i ; h α and h i are constant. together and analyzed according to the more restrictive scale: h p /h i .) The L -errors calculated at 0.3ssimulated time for the first group of grids — with h α varying between . cm and cm, h p /h i fixed at 4, and h i fixed at 1.1cm — are shown in Figure C.1c and have clear O ( h α ) convergence. The L -errors calculatedat 0.3s simulated time for the second group of grids — with h i varying between . cm and cm, h p /h i fixed at 4, and h α fixed at 3.3cm — are shown in Figure C.1d and have clear O ( h i ) convergence. Thisconvergence rate exceeds our expectations, but is likely an artifact of the regular Cartesian grid we use inthis particular simulation. And finally, the L -errors calculated at 1.28s simulated time for the final groupof grids — with h p /h i varying between 1 and 4, h i and h α fixed at 3.3cm — are shown in Figure C.1e andhave clear O ( h p /h i ) convergence.This validation of our method and implementation using the MMS shows that FV-MPM can reliablysolve the governing equations in (3.1) with at least first-order convergence in each discretization length scale.This low-order of convergence suggests that the method may be limited computationally for large complexflows; however, it does guarantee convergence for a wide range of problems which are extremely difficult toaddress with other numerical approaches (e.g. granular separation, solid-to-fluid and fluid-to-solid transition,long-duration simulations, etc.). Appendix D. Notes on Numerical Stability
In this section, we briefly discuss the stability of the numerical method described in the main text. Thereare two primary factors that affect this stability for the classes of engineering problems we are interested inaddressing: wave-propagation and dissipation. Assuming that the relevant features of a mixed flow solutionare properly resolved, we can express the limitations on our discrete time increment ( ∆ t ) that are necessaryto ensure stability when an explicit time integration approach is taken.39 ppendix D.1. Courant–Friedrichs–Lewy Condition The first condition is the Courant–Friedrichs–Lewy (CFL; see [87]) condition. As stated in the maintext, this condition generally requires that the characteristic lines associated with a hyperbolic system ofequations are appropriately resolved by the stencil of the numerical algorithm. In other words, the timeincrement must be small enough that the distance traveled by a physical wave (e.g. acoustic or elastic waves)does not exceed the relevant grid discretization length in a single time step. For the mixture problems weare generally interested in, this condition can be expressed as follows: ∆ t ≤ c h α (cid:13)(cid:13) a f + (cid:107) v f (cid:107) (cid:13)(cid:13) ∞ , and ∆ t ≤ c h i (cid:13)(cid:13) a s + (cid:107) v s (cid:107) (cid:13)(cid:13) ∞ , (D.1)where a f and a s are the acoustic wave speeds in the fluid and solid phase, respectively; c and c are O (1) constants; h i is the characteristic length of the finite element grid; and h α is the characteristic length of thefinite volume grid. In our simulations, we generally assume c , c (cid:46) . Appendix D.2. Effective Fluid Shear Stress
The second condition relates to the ability of our time integration procedure to resolve the rates ofdissipation in a physical mixture. One contributing factor to dissipation in our mixed flow problems is theviscous shear stress in the fluid, τ f . Although there many admissible forms of this stress, we will focus hereon those that obey the relationship in (2.10): τ f = η r (cid:0) ∇ v f + ∇ v (cid:62) f − div( v f ) I (cid:1) . When implemented numerically, as needed to solve (3.4), the shear stress is approximated along the finitevolume boundaries using either cell-centered or boundary-centered reconstructions. In either case, we cangenerally express the integrated shear stress along the boundary ∂ Ω ( β,γ ) (i.e. the cell boundary betweenvolumes Ω β and Ω γ ) as follows, (cid:90) ∂ Ω ( β,γ ) τ f ˆ n β → γ da ≈ η r dA ( β,γ ) ( C ( β,γ ) − C ( γ,β ) ) (cid:18) (cid:104) ¯ ρ f v f (cid:105) γ (cid:104) ¯ ρ f (cid:105) γ − (cid:104) ¯ ρ f v f (cid:105) β (cid:104) ¯ ρ f (cid:105) β (cid:19) + η r dA ( β,γ ) (cid:88) α ∈ S β / { γ } C ( β,α ) (cid:18) (cid:104) ¯ ρ f v f (cid:105) α (cid:104) ¯ ρ f (cid:105) α − (cid:104) ¯ ρ f v f (cid:105) β (cid:104) ¯ ρ f (cid:105) β (cid:19) + η r dA ( β,γ ) (cid:88) α ∈ S γ / { β } C ( γ,α ) (cid:18) (cid:104) ¯ ρ f v f (cid:105) α (cid:104) ¯ ρ f (cid:105) α − (cid:104) ¯ ρ f v f (cid:105) γ (cid:104) ¯ ρ f (cid:105) γ (cid:19) , (D.2)with ˆ n β → γ the oriented face normal pointing from Ω β to Ω γ ; dA ( β,γ ) the area of the boundary; C ( β,γ ) a gridspecific parameter used to calculate gradients between the centroids of Ω β and Ω γ (typically scales with h − α );and S β and S γ the set of neighbor indices for Ω β and Ω γ , respectively. Here ( C ( β,γ ) − C ( γ,β ) ) is non-zeroand positive.When substituted into (3.4), (D.2) represents one component of the momentum flux between the volumes Ω β and Ω γ . For positive ( C ( β,γ ) − C ( γ,β ) ) , this momentum flux will generally reduce the relative velocitiesbetween these neighboring cells; however, when an explicit time integrator (e.g. Forward Euler) is used witha large time increment ( ∆ t ), this diffusive flux can ‘over correct’ velocity differences and become unstable.To ensure (cid:104) ¯ ρ f v f (cid:105) β / (cid:104) ¯ ρ f (cid:105) β − (cid:104) ¯ ρ f v f (cid:105) γ / (cid:104) ¯ ρ f (cid:105) γ doesn’t explode numerically, we require that, ∆ t (cid:46) c h α η r (cid:2) (cid:104) ¯ ρ f (cid:105) − β + (cid:104) ¯ ρ f (cid:105) − γ (cid:3) , ∀ γ ∈ S β , ∀ β ∈ [1 , N v ] , (D.3)with c an O (2) constant and h α the characteristic length of the finite volume grid. In our simulations, wegenerally assume c (cid:46) . 40 ppendix D.3. Effective Fluid Heat Flux Another contributing factor to dissipation in mixed flow problems is the diffusion of temperature throughthe heat flux q f . Although there many admissible forms of this heat flux, we will focus here on those thatobey the relationship in (2.12): q f = − nk f ∇ ϑ f . When implemented numerically, as needed to solve (3.4), the heat flux is approximated along the finitevolume boundaries using either cell-centered or boundary-centered reconstructions of ∇ ϑ f . For the idealgas law considered in this work, we can express the reconstructed heat flux along the boundary ∂ Ω ( β,γ ) asfollows, (cid:90) ∂ Ω ( β,γ ) − q f · ˆ n β → γ da ≈ nk f c v dA ( β,γ ) ( C ( β,γ ) − C ( γ,β ) ) (cid:0) (cid:104) ε f (cid:105) γ − (cid:104) ε f (cid:105) β (cid:1) + nk f c v dA ( β,γ ) (cid:88) α ∈ S β / { γ } C ( β,α ) (cid:0) (cid:104) ε f (cid:105) α − (cid:104) ε f (cid:105) β (cid:1) + nk f c v dA ( β,γ ) (cid:88) α ∈ S γ / { β } C ( γ,α ) (cid:0) (cid:104) ε f (cid:105) α − (cid:104) ε f (cid:105) γ (cid:1) , (cid:104) ε f (cid:105) α ≈ (cid:104) ¯ ρ f (cid:105) α (cid:104) ¯ ρ f E f (cid:105) α − (cid:104) ¯ ρ f v f (cid:105) α (cid:104) ¯ ρ f (cid:105) α , ∀ α ∈ [1 , N v ] , (D.4)with c v the specific heat of the gas at constant volume and ˆ n β → γ , dA ( β,γ ) , C ( β,γ ) , and S β and S γ defined inthe previous section.When substituted into (3.4), (D.4) represents one component of the energy flux between volumes Ω β and Ω γ . For positive (cid:0) C ( β,γ ) − C ( γ,β ) (cid:1) , this energy flux will generally reduce the difference between the specificinternal energies of the neighboring cells; however, when an explicit time integrator (e.g. Forward Euler) isused with a large time increment, this diffusive term can ‘over correct’ the energy difference and becomeunstable. To ensure (cid:104) ε f (cid:105) β − (cid:104) ε f (cid:105) γ doesn’t explode numerically, we require that, ∆ t (cid:46) c c v h α nk f (cid:2) (cid:104) ¯ ρ f (cid:105) − β + (cid:104) ¯ ρ f (cid:105) − γ (cid:3) , ∀ γ ∈ S β , ∀ β ∈ [1 , N v ] , (D.5)with c and O (2) constant and h α the characteristic length of the finite volume grid. Appendix D.4. Inter-phase Drag Force
The final form of dissipation that can affect the stability of our method is the inter-phase drag, f d ,which tends to resist relative motion of the two mixed phases. An exact assessment of the stability conditionrelating to drag is beyond the scope of this Appendix; however, if we analyze the behavior of the approximatedrag force vectors in (B.6) and (B.7), f drag i ∗ = − K ∗ i (cid:18) V i v si − N v (cid:88) α =1 V α A iα (cid:104) ¯ ρ f v f (cid:105) α / (cid:104) ¯ ρ f (cid:105) α (cid:19) , F drag α ∗ = N n (cid:88) i =1 A iα K ∗ i (cid:0) v si − (cid:104) ¯ ρ f v f (cid:105) α / (cid:104) ¯ ρ f (cid:105) α (cid:1) · v si , with K ∗ i defined in (B.5), we can make a broad statements about when this term may become unstable andpropose an adjustment to enhance its stability. 41onsider the numerical expressions for momentum conservation in the solid and fluid phases from (3.16)and (3.17). If we limit our analysis to methods that uses the diagonal matrices [ M D ] and [ B D ] , then wehave: m i D s v s Dt = − K ∗ i (cid:18) V i v si − N v (cid:88) α =1 V α A iα (cid:104) ¯ ρ f v f (cid:105) α / (cid:104) ¯ ρ f (cid:105) α (cid:19) + f int i + f ext i + f buoy i + f τ i , (D.6)and, ddt (cid:32) (cid:104) ¯ ρ f (cid:105) α (cid:104) ¯ ρ f v f (cid:105) α (cid:104) ¯ ρ f E f (cid:105) α (cid:33) = N n (cid:88) i =1 A iα K ∗ i (cid:0) v si − (cid:104) ¯ ρ f v f (cid:105) α / (cid:104) ¯ ρ f (cid:105) α (cid:1) · (cid:32) v si (cid:33) + F int α + F ext α + F buoy α . (D.7)Integrating (D.6) and (D.7) over time will generally reduce the difference between the predicted solid phaseand fluid phase velocities; however, when an explicit time integrator (e.g. Forward Euler) is used with a largetime increment, these drag terms can ‘over correct’ the velocity difference and produce an unstable scheme.To ensure v s − v f doesn’t explode numerically, we require that, ∆ t ≤ K ∗ i (cid:2) ( φ i ρ s ) − − (cid:104) ¯ ρ f (cid:105) − α (cid:3) , ∀ ( i, α ) ∈ (cid:8) ( j, β ) | A jβ (cid:54) = 0 (cid:9) , (D.8)with φ i = 1 − n i .On the other hand, if (D.6) and (D.7) are integrated with an implicit time integrator (e.g. BackwardEuler), substantially larger time increments could be taken without incurring an ‘over correction’. Such atime integration approach generally requires the construction and inversion of system spanning matrices;however, in the special case of the inter-phase drag, this inversion can be performed independently of theother forcing terms. This suggests that a semi-implicit time integrator can be developed that allows forstable time integration without requiring a complicated matrix inversion.Here we propose such an integrator. (Note that a version of this semi-implicit drag is used in thealternative numerical algorithm described in the next section and in the simulations reported in sections 8.3,8.4, 8.5, and 8.6.) This approach uses explicit time integration to determine an intermediate mixture state( (cid:104) ¯ ρ f (cid:105) ∗ α , (cid:104) ¯ ρ f v f (cid:105) ∗ α , (cid:104) ¯ ρ f E f (cid:105) ∗ α , and v ∗ si ) between the time increments s and s + 1 : (cid:32) (cid:104) ¯ ρ f (cid:105) ∗ α (cid:104) ¯ ρ f v f (cid:105) ∗ α (cid:104) ¯ ρ f E f (cid:105) ∗ α (cid:33) := (cid:32) (cid:104) ¯ ρ f (cid:105) sα (cid:104) ¯ ρ f v f (cid:105) sα (cid:104) ¯ ρ f E f (cid:105) sα (cid:33) + ∆ t (cid:2) ( F int α ) s + ( F ext α ) s + ( F buoy α ) s (cid:3) ,m i v ∗ si := m i v ssi + ∆ t (cid:2) ( f int i ) s + ( f ext i ) s + ( f buoy i ) s + ( f τ i ) s (cid:3) , (D.9)with ( · ) s denoting an evaluation of the forcing vectors at the s th time increment. With these intermediatestates determined, the drag force vectors ( f drag i ∗ ) s +1 and ( F drag α ∗ ) s +1 can be determined as follows, ( f drag i ∗ ) s +1 = − ˜ K ∗ i (cid:18) V i v ∗ si − N v (cid:88) α =1 V α A iα (cid:104) ¯ ρ f v f (cid:105) ∗ α / (cid:104) ¯ ρ f (cid:105) ∗ α (cid:19) , (D.10) ( F drag α ∗ ) s +1 = N n (cid:88) i =1 A iα ˜ K ∗ i (cid:0) v ∗ si − (cid:104) ¯ ρ f v f (cid:105) ∗ α / (cid:104) ¯ ρ f (cid:105) ∗ α (cid:1) · v ∗ si , (D.11)with, ˜ K ∗ i ≈ ( K ∗ i ) s t ( K ∗ i ) s (cid:2) min α ∈V i ( (cid:104) ¯ ρ f (cid:105) α ) − + ( φ i ρ s ) − (cid:3) , V i = { α | A iα (cid:54) = 0 } . (D.12)This formulation is exact for locally uniform flows and ensures that the drag term remains stable regardlessof the chosen time increment. 42 ppendix E. Alternative Numerical Algorithm In this final section, we propose an alternative numerical algorithm for integrating the time-history ofthe solution coefficients in (3.23). (Note that this algorithm was used to generate the results describedin sections 8.3, 8.4, 8.5, and 8.6.) We use a first-order, Forward Euler time integration method (with theupdate-stress-last MPM approach) for the solid phase equations of motion, we use a higher-order integrationmethod (e.g. 4th-order Runge–Kutta) for the fluid phase equations, and we implement the semi-implicitdrag forces described in the previous section. The benefits of this approach are a larger numerical stencil(particularly for problems where the limiting time-scale is associated with the propagation of acoustic wavesin the fluid), increased accuracy in regions of the mixture without granular material, and the ability to uselarger time increments when the drag forces in the mixture are relatively large.Suppose the following parameters are known at t k : (i) the mapping matrices, [ S ] k , [ G ] k , [ A ] , [ B D ] , and [ M D ] k ; (ii) the solution coefficients, ¯ ρ ksp , ˜ σ kp , ¯ ξ kp , (cid:104) ¯ ρ f (cid:105) kα , (cid:104) ¯ ρ f v f (cid:105) kα , and (cid:104) ¯ ρ f E f (cid:105) kα ; (iii) the coefficients of thematerial point approximation of the velocity field, v ∗ ksp ; and (iv) the numerical representation of the materialpoint characteristic functions, U p ( x , t k ) , (e.g. v kp and x kp for basic MPM [28] and uGIMP [39]). We determinethe value of the solution coefficient at t k +1 according to the following steps:(1) Determine the mixture porosity coefficients, as in (5.3), using the diagonal matrix [ B D ] : B Dii (1 − n ki ) ρ s = N m (cid:88) p =1 S kip m p ∀ i ∈ [1 , N n ] . (E.1)(2) Determine the solid phase velocity coefficients, as in (5.6), using the diagonal matrix [ M D ] : M Dkii v ksi = N p (cid:88) p =1 S kip m p v ∗ ksp , ∀ i ∈ [1 , N n ] . (E.2)(3) Calculate the nodal force vectors, ( f int i ) k and ( f ext i ) k , as in (3.16).(4) Estimate ( f drag i ∗ ) k ; ( F drag α ∗ ) k ; and the full-step force vectors ( f buoy i ∗ ) k , ( F int α ) k , ( F ext α ) k , and ( F buoy α ) k :i) Construct an approximation of the fluid phase fields, ρ f , v f , E f , and ¯ ρ f , within each finite volumecell, as in (4.4).ii) Use the reconstructed fluid phase fields to calculate the nodal force vectors, ( f drag i ∗ ) k and ( f buoy i ∗ ) k ,as in (B.6) and (B.1), and the finite volume force vectors, ( F int α ) k , ( F ext α ) k , ( F drag α ∗ ) k , and ( F buoy α ) k , as in (3.17), (3.19), and (B.7).(5) Estimate the half-step force vectors ( f buoy i ∗ ) k , ( F int α ) k , ( F ext α ) k , and ( F buoy α ) k :i) Calculate intermediate fluid phase field coefficients from the equations of motion as: (cid:104) ¯ ρ f (cid:105) kα (cid:104) ¯ ρ f v f (cid:105) kα (cid:104) ¯ ρ f E f (cid:105) kα = (cid:104) ¯ ρ f (cid:105) kα (cid:104) ¯ ρ f v f (cid:105) kα (cid:104) ¯ ρ f E f (cid:105) kα + ∆ t (cid:2) ( F int α ) k +( F ext α ) k +( F buoy α ) k +( F drag α ∗ ) k (cid:3) , ∀ α ∈ [1 , N v ] . (E.3)ii) Construct an intermediate approximation of the fluid phase fields, ρ f , v f , E f , and ¯ ρ f , withineach finite volume cell, using {(cid:104) ¯ ρ f (cid:105) kα } , {(cid:104) ¯ ρ f v f (cid:105) kα } , and {(cid:104) ¯ ρ f E f (cid:105) kα } .iii) Use the reconstructed fluid phase fields to calculate the nodal force vector ( f buoy i ∗ ) k , as in (B.1),and the finite volume force vectors, ( F int α ) k , ( F ext α ) k , and ( F buoy α ) k , as in (3.17) and (3.19).(6) Estimate the half-step force vectors ( f buoy i ∗ ) k , ( F int α ) k , ( F ext α ) k , and ( F buoy α ) k :43) Calculate intermediate fluid phase field coefficients from the equations of motion as: (cid:104) ¯ ρ f (cid:105) kα (cid:104) ¯ ρ f v f (cid:105) kα (cid:104) ¯ ρ f E f (cid:105) kα = (cid:104) ¯ ρ f (cid:105) kα (cid:104) ¯ ρ f v f (cid:105) kα (cid:104) ¯ ρ f E f (cid:105) kα + ∆ t (cid:2) ( F int α ) k +( F ext α ) k +( F buoy α ) k +( F drag α ∗ ) k (cid:3) , ∀ α ∈ [1 , N v ] . (E.4)ii) Construct the second intermediate approximation of the fluid phase fields, ρ f , v f , E f , and ¯ ρ f ,within each finite volume cell, using {(cid:104) ¯ ρ f (cid:105) kα } , {(cid:104) ¯ ρ f v f (cid:105) kα } , and {(cid:104) ¯ ρ f E f (cid:105) kα } .iii) Use the reconstructed fluid phase fields to calculate the nodal force vector ( f buoy i ∗ ) k , as in (B.1),and the finite volume force vectors, ( F int α ) k , ( F ext α ) k , and ( F buoy α ) k , as in (3.17) and (3.19).(7) Estimate the full-step force vectors ( f buoy i ∗ ) k , ( F int α ) k , ( F ext α ) k , and ( F buoy α ) k :i) Calculate the third intermediate fluid phase field coefficients from the equations of motion as: (cid:104) ¯ ρ f (cid:105) kα (cid:104) ¯ ρ f v f (cid:105) kα (cid:104) ¯ ρ f E f (cid:105) kα = (cid:104) ¯ ρ f (cid:105) kα (cid:104) ¯ ρ f v f (cid:105) kα (cid:104) ¯ ρ f E f (cid:105) kα +∆ t (cid:2) ( F int α ) k +( F ext α ) k +( F buoy α ) k +( F drag α ∗ ) k (cid:3) , ∀ α ∈ [1 , N v ] . (E.5)ii) Construct the third intermediate approximation of the fluid phase fields, ρ f , v f , E f , and ¯ ρ f ,within each finite volume cell, using {(cid:104) ¯ ρ f (cid:105) kα } , {(cid:104) ¯ ρ f v f (cid:105) kα } , and {(cid:104) ¯ ρ f E f (cid:105) kα } .iii) Use the reconstructed fluid phase fields to calculate the nodal force vector ( f buoy i ∗ ) k , as in (B.1),and the finite volume force vectors, ( F int α ) k , ( F ext α ) k , and ( F buoy α ) k , as in (3.17) and (3.19).(8) Estimate the semi-implicit force vectors ( f drag i ∗ ) k +1 and ( F drag α ∗ ) k +1 :i) Calculate the intermediate fluid field coefficients for the semi-implicit calculation in (D.9): (cid:32) (cid:104) ¯ ρ f (cid:105) ∗ α (cid:104) ¯ ρ f v f (cid:105) ∗ α (cid:104) ¯ ρ f E f (cid:105) ∗ α (cid:33) = (cid:104) ¯ ρ f (cid:105) kα (cid:104) ¯ ρ f v f (cid:105) kα (cid:104) ¯ ρ f E f (cid:105) kα + ∆ t (cid:2) ( F int α ) k + ( F ext α ) k + ( F buoy α ) k (cid:3) + ∆ t (cid:2) ( F int α ) k + ( F ext α ) k + ( F buoy α ) k (cid:3) + ∆ t (cid:2) ( F int α ) k + ( F ext α ) k + ( F buoy α ) k (cid:3) + ∆ t (cid:2) ( F int α ) k + ( F ext α ) k + ( F buoy α ) k (cid:3) , ∀ α ∈ [1 , N v ] . (E.6)ii) Calculate the intermediate solid velocity coefficients for the semi-implicit calculation in (D.9): M Dkii v ∗ si = M Dkii v ksi + ∆ t (cid:2) ( f int i ) k + ( f ext i ) k + f τ i (cid:3) + ∆ t (cid:2) ( f buoy i ∗ ) k + 2( f buoy i ∗ ) k + 2( f buoy i ∗ ) k + ( f buoy i ∗ ) k (cid:3) , ∀ i ∈ [1 , N n ] . (E.7)iii) Use these intermediate coefficients to calculate the nodal force vector ( f drag i ∗ ) k +1 , as in (D.10),and the finite volume force vector ( F drag α ∗ ) k +1 , as in (D.11).449) Now calculate the updated fluid phase field coefficients from the equations of motion as: (cid:104) ¯ ρ f (cid:105) k +1 α (cid:104) ¯ ρ f v f (cid:105) k +1 α (cid:104) ¯ ρ f E f (cid:105) k +! α = (cid:104) ¯ ρ f (cid:105) kα (cid:104) ¯ ρ f v f (cid:105) kα (cid:104) ¯ ρ f E f (cid:105) kα + ∆ t (cid:2) ( F int α ) k + ( F ext α ) k + ( F buoy α ) k + ( F drag α ∗ ) k +1 (cid:3) + ∆ t (cid:2) ( F int α ) k + ( F ext α ) k + ( F buoy α ) k + ( F drag α ∗ ) k +1 (cid:3) + ∆ t (cid:2) ( F int α ) k + ( F ext α ) k + ( F buoy α ) k + ( F drag α ∗ ) k +1 (cid:3) + ∆ t (cid:2) ( F int α ) k + ( F ext α ) k + ( F buoy α ) k + ( F drag α ∗ ) k +1 (cid:3) , ∀ α ∈ [1 , N v ] . (E.8)(10) Determine the solid phase acceleration coefficients associated with the finite element grid nodes fromthe equations of motion and appropriate boundary conditions: M Dkii a ksi = ( f int i ) k + ( f ext i ) k + ( f drag i ∗ ) k +1 + f τ i + 16 (cid:2) ( f buoy i ∗ ) k + 2( f buoy i ∗ ) k + 2( f buoy i ∗ ) k + ( f buoy i ∗ ) k (cid:3) , ∀ i ∈ [1 , N n ] . (E.9)(11) Use an explicit time integrator to obtain the updated solid phase velocity coefficients associated withthe finite element grid nodes; treat the grid as though it were Lagrangian: v k +1 (cid:48) si = v ksi + ∆ t a ksi , ∀ i ∈ [1 , N n ] . (E.10)(12) Update the solid phase material state vector, ¯ ξ p , and effective granular stress, ˜ σ p , according to therelevant constitutive update procedure. For stability, L p , the average material point velocity gradientis often used: L k +1 p = N n (cid:88) i =1 v k +1 (cid:48) si ⊗ G kip , ˜ σ k +1 p = T (cid:0) L k +1 p , ¯ ξ k +1 p (cid:1) , ∀ p ∈ [1 , N m ] . (E.11)(13) Map the solid phase nodal accelerations and velocities to the material points to update their velocityapproximations, positions, and densities: v ∗ k +1 sp = v ∗ ksp + ∆ t N n (cid:88) i =1 S kip a ksi , ∀ p ∈ [1 , N m ] , x k +1 p = x kp + ∆ t N n (cid:88) i =1 S kip v k +1 (cid:48) si , ∀ p ∈ [1 , N m ] ,v k +1 p = v kp (cid:18) t N n (cid:88) i =1 G kip · v k +1 (cid:48) si (cid:19) , ∀ p ∈ [1 , N m ] , ¯ ρ k +1 sp = m p /v k +1 p , ∀ p ∈ [1 , N m ] . (E.12)(14) Update the diagonal mass matrix, [ M D ] k +1 , and mapping matrices, [ S ] k +1 and [ G ] k +1 , according totheir definitions in (3.13).(15) Increment k to k + 1 ; go to (1).This algorithm is consistent with the governing equations in (3.1) and has similar stability conditions tothe algorithm in the main text; however, the primary advantage of this approach is improved stability in thefluid phase of the simulation and increased accuracy in regions of the mixed flow that are without granularmaterial. 45 ppendix F. Artificial Viscosity for Strong Shocks The use of artificial (i.e. numerical) viscosity is common in gas dynamics simulations for problems in-volving strong shocks (see [107]). For the rocket exhaust impingement simulation in the main text, we usea simple shock-capturing scheme based on the works of [105, 106, 108]. In these approaches, the numericalequations for the fluid phase of the mixture is augmented as follows: ddt (cid:32) (cid:104) ¯ ρ f (cid:105) α (cid:104) ¯ ρ f v f (cid:105) α (cid:104) ¯ ρ f E f (cid:105) α (cid:33) = F int α + F ext α + F drag α + F buoy α + F vα . (F.1)Here F vα represents the additional flux arising from the artificial viscosity and takes the following form: F vα = (cid:90) ∂ Ω α (cid:15) v h ( u + − u − ) da, (F.2)with u + and u − the two reconstructed fluid state vectors on each side of the boundary, ∂ Ω α ; (cid:15) v the scalarartificial viscosity; and h the local grid length scale. (Note that (cid:15) v only appears in the fraction (cid:15) v /h and isdefined below.) In our implementation we use the following shock capturing form of (cid:15) v : (cid:15) v h = λ max min (cid:18) ( hλ ) β a min , (cid:19) ,hλ = (cid:40) | a + f − a − f | − ( v + f − v − f ) · ˆ n , if ( v + f − v − f ) · ˆ n < , else , (F.3)with β = 0 . ; ˆ n the oriented face normal; v + f and v − f the reconstructed fluid velocities on each side of theboundary; a + f and a − f the local acoustic wave speeds associated with these reconstructions; a min the smallerof these two values; and λ max = max (cid:0) a + + (cid:107) v + f (cid:107) , a − + (cid:107) v − f (cid:107) (cid:1) . (F.4)46 eferences [1] Raffi M Turian and Tran-Fu Yuan. Flow of slurries in pipelines. AIChE Journal , 23(3):232–243, 1977.[2] Mickael Pailha and Olivier Pouliquen. A two-phase flow description of the initiation of underwatergranular avalanches.
Journal of Fluid Mechanics , 633:115–135, 2009.[3] Jeffrey D Keller, Glen R Whitehouse, Daniel A Wachspress, Milton E Teske, and Todd R Quacken-bush. A physics-based model of rotorcraft brownout for flight simulation applications. In
ANNUALFORUM PROCEEDINGS-AMERICAN HELICOPTER SOCIETY , volume 62, page 1098. AMERI-CAN HELICOPTER SOCIETY, INC, 2006.[4] Shashank Agarwal, Carmine Senatore, Tingnan Zhang, Mark Kingsbury, Karl Iagnemma, Daniel IGoldman, and Ken Kamrin. Modeling of the interaction of rigid wheels with dry granular media.
Journal of Terramechanics , 85:1–14, 2019.[5] Albert Einstein. Calculation of the viscosity-coefficient of a liquid in which a large number of smallspheres are suspended in irregular distribution.
Ann. Phys. Leipzig , 19:286–306, 1906.[6] GK Batchelor and JT Green. The determination of the bulk stress in a suspension of spherical particlesto order c 2.
Journal of Fluid Mechanics , 56(3):401–427, 1972.[7] JS Chong, EB Christiansen, and AD Baer. Rheology of concentrated suspensions.
Journal of appliedpolymer science , 15(8):2007–2021, 1971.[8] Philip Crosbie Carman. Fluid flow through granular beds.
Trans. Inst. Chem. Eng. , 15:150–166, 1937.[9] Benjamin K Cook, David R Noble, and John R Williams. A direct simulation method for particle-fluidsystems.
Engineering Computations , 21(2/3/4):151–168, 2004.[10] ZY Zhou, SB Kuang, KW Chu, and AB Yu. Discrete particle simulation of particle-fluid flow: modelformulations and their applicability.
Journal of Fluid Mechanics , 661:482, 2010.[11] Ryohei Seto, Romain Mari, Jeffrey F Morris, and Morton M Denn. Discontinuous shear thickening offrictional hard-sphere suspensions.
Physical review letters , 111(21):218301, 2013.[12] Romain Mari, Ryohei Seto, Jeffrey F Morris, and Morton M Denn. Discontinuous shear thickening inbrownian suspensions by dynamic simulation.
Proceedings of the National Academy of Sciences , 112(50):15326–15330, 2015.[13] Lhassan Amarsid, J-Y Delenne, Patrick Mutabaruka, Yann Monerie, Frédéric Perales, and FarhangRadjai. Viscoinertial regime of immersed granular flows.
Physical Review E , 96(1):012901, 2017.[14] Jia Mao, Lanhao Zhao, Yingtang Di, Xunnan Liu, and Weiya Xu. A resolved cfd-dem approach for thesimulation of landslides and impulse waves.
Computer Methods in Applied Mechanics and Engineering ,359:112750, 2020. ISSN 0045-7825. doi: https://doi.org/10.1016/j.cma.2019.112750.[15] Ning Guo and Jidong Zhao. Parallel hierarchical multiscale modelling of hydro-mechanical problemsfor saturated granular soils.
Computer Methods in Applied Mechanics and Engineering , 305:37 – 61,2016. ISSN 0045-7825. doi: https://doi.org/10.1016/j.cma.2016.03.004.[16] Krishna Kumar, J-Y Delenne, and Kenichi Soga. Mechanics of granular column collapse in fluid atvarying slope angles.
Journal of Hydrodynamics , 29(4):529–541, 2017.[17] Kenichi Soga, E Alonso, A Yerro, K Kumar, and S Bandara. Trends in large-deformation analysis oflandslide mass movements with particular emphasis on the material point method.
Géotechnique , 66(3):248–273, 2015. 4718] Francesca Ceccato and Paolo Simonini. Granular flow impact forces on protection structures: Mpmnumerical simulations with different constitutive models.
Procedia Engineering , 158:164–169, 2016.[19] Elliot James Fern and Kenichi Soga. The role of constitutive models in mpm simulations of granularcolumn collapses.
Acta Geotechnica , 11(3):659–678, 2016.[20] Samila Bandara and Kenichi Soga. Coupling of soil deformation and pore fluid flow using materialpoint method.
Computers and geotechnics , 63:199–214, 2015.[21] Aaron S Baumgarten and Ken Kamrin. A general fluid–sediment mixture model and constitutivetheory validated in many flow regimes.
Journal of Fluid Mechanics , 861:721–764, 2019.[22] Olivier Coussy.
Poromechanics . John Wiley & Sons, 2004.[23] DS Drumheller. On theories for reacting immiscible mixtures.
International journal of engineeringscience , 38(3):347–382, 2000.[24] Roy Jackson.
The dynamics of fluidized particles . Cambridge University Press, 2000.[25] Václav Klika. A guide through available mixture theories for applications.
Critical Reviews in SolidState and Materials Sciences , 39(2):154–174, 2014.[26] Pascale Aussillous, Julien Chauchat, Mickael Pailha, Marc Médale, and Elisabeth Guazzelli. Investi-gation of the mobile granular layer in bedload transport by laminar shearing flows.
Journal of FluidMechanics , 736:594–615, 2013.[27] Ryan S Mieras, Jack A Puleo, Dylan Anderson, Tian-Jian Hsu, Daniel T Cox, and Joseph Calantoni.Relative contributions of bed load and suspended load to sediment transport under skewed-asymmetricwaves on a sandbar crest.
Journal of Geophysical Research: Oceans , 124(2):1294–1321, 2019.[28] Deborah Sulsky, Zhen Chen, and Howard L Schreyer. A particle method for history-dependent mate-rials.
Computer methods in applied mechanics and engineering , 118(1-2):179–196, 1994.[29] Keita Abe, Kenichi Soga, and Samila Bandara. Material point method for coupled hydromechanicalproblems.
Journal of Geotechnical and Geoenvironmental Engineering , 140(3):04013033, 2013.[30] Dongfang Liang, Xuanyu Zhao, and Mario Martinelli. Mpm simulations of the interaction betweenwater jet and soil bed.
Procedia Engineering , 175:242–249, 2017.[31] Irene Redaelli, Francesca Ceccato, Claudio Di Prisco, and Paolo Simonini. Solid-fluid transition ingranular flows: Mpm simulations with a new constitutive approach.
Procedia Engineering , 175:80–85,2017.[32] Francesca Ceccato, Lars Beuth, Pieter A Vermeer, and Paolo Simonini. Two-phase material pointmethod applied to the study of cone penetration.
Computers and Geotechnics , 80:440–452, 2016.[33] Aaron S. Baumgarten and Ken Kamrin. A general constitutive model for dense, fine-particle sus-pensions validated in many geometries.
Proceedings of the National Academy of Sciences , 116(42):20828–20836, 2019. ISSN 0027-8424. doi: 10.1073/pnas.1908065116. URL .[34] Wen-Chia Yang, Pedro Arduino, Gregory R Miller, and Peter Mackenzie-Helnwein. Smoothing algo-rithm for stabilization of the material point method for fluid–solid interaction problems.
ComputerMethods in Applied Mechanics and Engineering , 342:177–199, 2018.[35] SG Bardenhagen. Energy conservation error in the material point method for solid mechanics.
Journalof Computational Physics , 180(1):383–403, 2002.4836] M Steffen, PC Wallstedt, JE Guilkey, RM Kirby, and M Berzins. Examination and analysis of imple-mentation choices within the material point method (mpm).
Computer Modeling in Engineering andSciences , 32(2):107–127, 2008.[37] Michael Steffen, Robert M Kirby, and Martin Berzins. Decoupling and balancing of space and timeerrors in the material point method (mpm).
International journal for numerical methods in engineering ,82(10):1207–1243, 2010.[38] Deborah Sulsky and Ming Gong. Improving the material-point method. In
Innovative NumericalApproaches for Multi-Field and Multi-Scale Problems , pages 217–240. Springer, 2016.[39] Alireza Sadeghirad, Rebecca M Brannon, and Jeff Burghardt. A convected particle domain interpo-lation technique to extend applicability of the material point method for problems involving massivedeformations.
International Journal for numerical methods in Engineering , 86(12):1435–1456, 2011.[40] Yonggang Zheng, Fei Gao, Hongwu Zhang, and Mengkai Lu. Improved convected particle domaininterpolation method for coupled dynamic analysis of fully saturated porous media involving largedeformation.
Computer Methods in Applied Mechanics and Engineering , 257:150–163, 2013.[41] A Sadeghirad, Rebecca M Brannon, and JE Guilkey. Second-order convected particle domain interpo-lation (cpdi2) with enrichment for weak discontinuities at material interfaces.
International Journalfor numerical methods in Engineering , 95(11):928–952, 2013.[42] Yonghao Yue, Breannan Smith, Christopher Batty, Changxi Zheng, and Eitan Grinspun. Continuumfoam: A material point method for shear-dependent flows.
ACM Transactions on Graphics (TOG) ,34(5):160, 2015.[43] Vinh Phu Nguyen, Chi Thanh Nguyen, Timon Rabczuk, and Sundararajan Natarajan. On a family ofconvected particle domain interpolations in the material point method.
Finite Elements in Analysisand Design , 126:50–64, 2017.[44] Shyamini Kularathna and Kenichi Soga. Implicit formulation of material point method for analysisof incompressible materials.
Computer Methods in Applied Mechanics and Engineering , 313:673–686,2017.[45] Georgios Moutsanidis, Christopher C. Long, and Yuri Bazilevs. Iga-mpm: The isogeometric materialpoint method.
Computer Methods in Applied Mechanics and Engineering , 372:113346, 2020. ISSN0045-7825. doi: https://doi.org/10.1016/j.cma.2020.113346.[46] L GEORGE DAVID and M RICHARD. A two-phase debris-flow model that includes coupled evolutionof volume fractions, granular dilatancy, and pore-fluid pressure.
Italian journal of engineering geologyand Environment , 43:415–424, 2011.[47] Weiming Wu and Sam SY Wang. One-dimensional explicit finite-volume model for sediment transport.
Journal of Hydraulic Research , 46(1):87–98, 2008.[48] Mingliang Zhang and WM Wu. A two dimensional hydrodynamic and sediment transport model fordam break based on finite volume method with quadtree grid.
Applied Ocean Research , 33(4):297–308,2011.[49] Sebastian Geiger, Stephen Roberts, Stephan K Matthäi, Christopher Zoppou, and Adrian Burri. Com-bining finite element and finite volume methods for efficient multiphase flow simulations in highlyheterogeneous and structurally complex geologic media.
Geofluids , 4(4):284–299, 2004.[50] I. Rees, I. Masters, A.G. Malan, and R.W. Lewis. An edge-based finite volume scheme for saturated-unsaturated groundwater flow.
Computer Methods in Applied Mechanics and Engineering , 193(42):4741 – 4759, 2004. ISSN 0045-7825. doi: https://doi.org/10.1016/j.cma.2004.04.003.4951] Cyril W Hirt and Billy D Nichols. Volume of fluid (vof) method for the dynamics of free boundaries.
Journal of computational physics , 39(1):201–225, 1981.[52] James A Sethian and Peter Smereka. Level set methods for fluid interfaces.
Annual review of fluidmechanics , 35(1):341–372, 2003.[53] HS Udaykumar, R Mittal, P Rampunggoon, and A Khanna. A sharp interface cartesian grid method forsimulating flows with complex moving boundaries.
Journal of computational physics , 174(1):345–380,2001.[54] Lin Sun, Sanjay R Mathur, and Jayathi Y Murthy. An unstructured finite-volume method for incom-pressible flows with complex immersed boundaries.
Numerical Heat Transfer, Part B: Fundamentals ,58(4):217–241, 2010.[55] Philip L Roe. Approximate riemann solvers, parameter vectors, and difference schemes.
Journal ofcomputational physics , 43(2):357–372, 1981.[56] Anthony Jameson and D Mavriplis. Finite volume solution of the two-dimensional euler equations ona regular triangular mesh.
AIAA journal , 24(4):611–618, 1986.[57] Timothy Barth and Dennis Jespersen. The design and application of upwind schemes on unstructuredmeshes. In , page 366, 1989.[58] JY Trepanier, M Reggio, H Zhang, and R Camarero. A finite-volume method for the euler equationson arbitrary lagrangian-eulerian grids.
Computers & fluids , 20(4):399–409, 1991.[59] Aaron S Baumgarten. A coupled, two-phase fluid-sediment material model and mixture theory im-plemented using the material point method. Master’s thesis, Massachusetts Institute of Technology,2018.[60] Georgios Moutsanidis, David Kamensky, JS Chen, and Yuri Bazilevs. Hyperbolic phase field modelingof brittle fracture: Part ii—immersed iga–rkpm coupling for air-blast–structure interaction.
Journalof the Mechanics and Physics of Solids , 121:114–132, 2018.[61] Anvar Gilmanov and Sumanta Acharya. A hybrid immersed boundary and material point methodfor simulating 3d fluid–structure interaction problems.
International journal for numerical methods influids , 56(12):2151–2177, 2008.[62] ZP Chen, XM Qiu, X Zhang, and YP Lian. Improved coupling of finite element method with ma-terial point method based on a particle-to-surface contact algorithm.
Computer Methods in AppliedMechanics and Engineering , 293:1–19, 2015.[63] Yu Zhao and Robert H Davis. Interaction of two touching spheres in a viscous fluid.
Chemicalengineering science , 57(11):1997–2006, 2002.[64] K_H Roscoe and JB Burland. On the generalized stress-strain behaviour of wet clay.
EngineeringPlasticity , pages 535–609, 1968.[65] Karl Terzaghi.
Theoretical Soil Mechanics . John Wiley and Sons Inc., New York, New York, 1943.[66] Pierre Jop, Yoël Forterre, and Olivier Pouliquen. A constitutive law for dense granular flows.
Nature ,441(7094):727, 2006.[67] Sachith Dunatunga and Ken Kamrin. Continuum modelling and simulation of granular flows throughtheir many phases.
Journal of Fluid Mechanics , 779:483–513, 2015.[68] François Boyer, Élisabeth Guazzelli, and Olivier Pouliquen. Unifying suspension and granular rheology.
Physical Review Letters , 107(18):188301, 2011.5069] James T Jenkins and Stuart B Savage. A theory for the rapid flow of identical, smooth, nearly elastic,spherical particles.
Journal of fluid mechanics , 130:187–202, 1983.[70] Nikolai V Brilliantov and Thorsten Pöschel.
Kinetic theory of granular gases . Oxford University Press,2010.[71] Martin Anton van der Hoef, R Beetstra, and JAM Kuipers. Lattice-boltzmann simulations of low-reynolds-number flow past mono-and bidisperse arrays of spheres: results for the permeability anddrag force.
Journal of fluid mechanics , 528:233–254, 2005.[72] R Beetstra, Martin Anton van der Hoef, and JAM Kuipers. Drag force of intermediate reynolds numberflow past mono-and bidisperse arrays of spheres.
AIChE journal , 53(2):489–501, 2007.[73] Jules Étienne Juvénal Dupuit.
Études théoriques et pratiques sur le mouvement des eaux dans lescanaux découverts et à travers les terrains perméables: avec des considérations relatives au régime desgrandes eaux, au débouché à leur donner, et à la marche des alluvions dans les rivières à fond mobile .Dunod, 1863.[74] Jonathan J Stickel and Robert L Powell. Fluid mechanics and rheology of dense suspensions.
Annu.Rev. Fluid Mech. , 37:129–149, 2005.[75] GK Batchelor. The stress system in a suspension of force-free particles.
Journal of fluid mechanics ,41(3):545–570, 1970.[76] Krzysztof Wilmanski. Tortuosity and objective relative accelerations in the theory of porous materials.
Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences , 461(2057):1533–1561, 2005.[77] W Kosiński, J Kubik, and K Hutter. On the added mass effect for porous media.
Archives of Mechanics ,54(5-6):479–495, 2002.[78] Scott G Bardenhagen and Edward M Kober. The generalized interpolation material point method.
Computer Modeling in Engineering and Sciences , 5(6):477–496, 2004.[79] Michael Steffen, Robert M Kirby, and Martin Berzins. Analysis and reduction of quadrature errors inthe material point method (mpm).
International journal for numerical methods in engineering , 76(6):922–948, 2008.[80] Ami Harten. High resolution schemes for hyperbolic conservation laws.
Journal of computationalphysics , 135(2):260–278, 1997.[81] Jeremiah U Brackbill and Hans M Ruppel. Flip: A method for adaptively zoned, particle-in-cellcalculations of fluid flows in two dimensions.
Journal of Computational physics , 65(2):314–343, 1986.[82] Jeremiah U Brackbill, Douglas B Kothe, and Hans M Ruppel. Flip: a low-dissipation, particle-in-cellmethod for fluid flow.
Computer Physics Communications , 48(1):25–38, 1988.[83] Chenfanfu Jiang, Craig Schroeder, Andrew Selle, Joseph Teran, and Alexey Stomakhin. The affineparticle-in-cell method.
ACM Transactions on Graphics (TOG) , 34(4):51, 2015.[84] Chad C Hammerquist and John A Nairn. A new method for material point method particle updatesthat reduces noise and enhances stability.
Computer Methods in Applied Mechanics and Engineering ,318:724–738, 2017.[85] Thomas JR Hughes.
The Finite Element Method: Linear Static and Dynamic Finite Element Analysis .Prentice-Hall: Englewood Cliffs, 1987. 5186] Olivier Buzzi, Dorival M Pedroso, Anna Giacomini, et al. Caveats on the implementation of thegeneralized material point method.
Computer Modeling in Engineering and Sciences , 1(1):1–21, 2008.[87] Richard Courant, Kurt Friedrichs, and Hans Lewy. Über die partiellen differenzengleichungen dermathematischen physik.
Mathematische annalen , 100(1):32–74, 1928.[88] Henry Darcy.
Les fontaines publiques de la ville de Dijon . Victor Dalmont, 1856.[89] Karl Terzaghi, Ralph B Peck, and Gholamreza Mesri.
Soil mechanics in engineering practice . JohnWiley & Sons, 1996.[90] Loïc Rondon, Olivier Pouliquen, and Pascale Aussillous. Granular collapse in a fluid: role of the initialvolume fraction.
Physics of Fluids , 23(7):073301, 2011.[91] Christopher E Gritton.
Ringing instabilities in particle methods . PhD thesis, School of Computing,University of Utah, 2014.[92] Chris Gritton and Martin Berzins. Improving accuracy in the mpm method using a null space filter.
Computational Particle Mechanics , 4(1):131–142, 2017.[93] Fan Zhang, Xiong Zhang, and Yan Liu. An augmented incompressible material point method formodeling liquid sloshing problems.
International Journal of Mechanics and Materials in Design , 14(1):141–155, 2018.[94] Joseph Smagorinsky. General circulation experiments with the primitive equations: I. the basic exper-iment.
Monthly weather review , 91(3):99–164, 1963.[95] David C Wilcox et al.
Turbulence modeling for CFD , volume 2. La Canada, CA: DCM industries,1998.[96] Veit Schwämmle and Hans J Herrmann. Solitary wave behaviour of sand dunes.
Nature , 426(6967):619–620, 2003.[97] Gerd Sauermann, Pierre Rognon, A Poliakov, and Hans J Herrmann. The shape of the barchan dunesof southern morocco.
Geomorphology , 36(1-2):47–62, 2000.[98] Aaron B Morris, David B Goldstein, Philip L Varghese, and Laurence M Trafton. Approach formodeling rocket plume impingement and dust dispersal on the moon.
Journal of Spacecraft andRockets , 52(2):362–374, 2015.[99] Omid Ejtehadi and RS Myong. A modal discontinuous galerkin method for simulating dusty andgranular gas flows in thermal non-equilibrium in the eulerian framework.
Journal of ComputationalPhysics , page 109410, 2020.[100] PT Metzger, John E Lane, CD Immer, JN Gamsky, W Hauslein, Xiaoyi Li, RC Latta, III, andCM Donaue.
Scaling of erosion rate in subsonic jet experiments and Apollo lunar module landings ,pages 191–207. 2010.[101] Christopher Immer, Philip Metzger, Paul E Hintze, Andrew Nick, and Ryan Horan. Apollo 12 lunarmodule exhaust plume impingement on lunar surveyor iii.
Icarus , 211(2):1089–1102, 2011.[102] KL Farrow, JA German, and TM Gernstein. Operating characteristics of the apollo lm descent enginewith helium ingestion and propellant depletion in a simulated space environment. Technical report,ARNOLD ENGINEERING DEVELOPMENT CENTER ARNOLD AFB TN, 1967.[103] Jack M Cherne.
Mechanical Design of the Lunar Module Descent Engine . TRW Systems, 1967.52104] Howard A Perko, John D Nelson, and Jacklyn R Green. Mars soil mechanical properties and suitabilityof mars soil simulants.
Journal of Aerospace Engineering , 19(3):169–176, 2006.[105] Peter McCorquodale and Phillip Colella. A high-order finite-volume method for conservation laws onlocally refined grids.
Communications in Applied Mathematics and Computational Science , 6(1):1–25,2011.[106] Maurizio Pandolfi and Domenic D’Ambrosio. Numerical instabilities in upwind methods: analysis andcures for the “carbuncle” phenomenon.
Journal of Computational Physics , 166(2):271–301, 2001.[107] Mark L Wilkins. Use of artificial viscosity in multidimensional fluid dynamic calculations.
Journal ofcomputational physics , 36(3):281–303, 1980.[108] Garrett E Barter and David L Darmofal. Shock capturing with pde-based artificial viscosity for dgfem:Part i. formulation.
Journal of Computational Physics , 229(5):1810–1827, 2010.[109] Qiong Zhang and Ken Kamrin. Microscopic description of the granular fluidity field in nonlocal flowmodeling.
Physical review letters , 118(5):058001, 2017.[110] Seongmin Kim and Ken Kamrin. Power-law scaling in granular rheology across flow geometries.
Phys.Rev. Lett. , 125:088002, Aug 2020. doi: 10.1103/PhysRevLett.125.088002. URL https://link.aps.org/doi/10.1103/PhysRevLett.125.088002 .[111] P Mackenzie-Helnwein, P Arduino, W Shin, JA Moore, and GR Miller. Modeling strategies for multi-phase drag interactions using the material point method.
International journal for numerical methodsin engineering , 83(3):295–322, 2010.[112] Philip C Wallstedt and JE Guilkey. An evaluation of explicit time integration schemes for use withthe generalized interpolation material point method.