Coupling Particle Simulations: Challenges, Strategies, and the ON-DEM Vision
Full text
Coupling Particle Simulations Challenges, Strategies, and the ON-DEM Vision Thornton, Weinhart, Scrase, Schneider, Post, Barreto, et al 6th preCiCE Workshop, 11 September 2025
Open-source coupling with 1 / 30
Grains in industry 2 / 30
Grains in nature Tephra expulsion Flow through contraction Segregation Elongated runoff Fingering instability 3 / 30
The force model •Discrete particle model governed by Newtonian mechanics: mi d2xi dt2= fi •Contact forces and body forces: fi=X j fij + bi, •Contact force model: fij =fn ijn +ft ij t, fn ij =kδij+γvn ij, ft ij =−min(µfn ij, ktδt ij+γtvt ij) (Luding, 2008, Enviro. and Civil. Eng.) 4 / 30
The force model •Discrete particle model governed by Newtonian mechanics: mi d2xi dt2= fi •Contact forces and body forces: fi=X j fij + bi, •Contact force model: fij =fn ijn +ft ij t, fn ij =kδij+γvn ij, ft ij =−min(µfn ij, ktδt ij+γtvt ij) (Luding, 2008, Enviro. and Civil. Eng.) 4 / 30
The force model •Discrete particle model governed by Newtonian mechanics: mi d2xi dt2= fi •Contact forces and body forces: fi=X j fij + bi, •Contact force model: fij =fn ijn +ft ij t, fn ij =kδij+γvn ij, ft ij =−min(µfn ij, ktδt ij+γtvt ij) (Luding, 2008, Enviro. and Civil. Eng.) 4 / 30
Particle simulations with MercuryDPM •Fast Fast contact detection for polydisperse simulations •Flexible Complex walls and boundary conditions •Accurate Coarse-graining technique to extract continuum fields •Open-source Available at www.mercurydpm.org •Coupled/interfaced with , , , , GrainLearning y [m] x [m] 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0 mm 1 mm 2 mm 3 mm 4 mm 5 / 30
Fast computations, optimal contact detection Traditional MercuryDPM Speed-Up Linked Cell Hierarchical Grid for high polydispersity #checks: 1000 #checks: 89 up to 200x quicker! 6 / 30
Coupling methods Concurrent multi-scale modeling of particle-continuum interactions 3 (a) (b) (c) Fig. 1: (a) 2D illustration of the discrete particle and finite element domains WDE and WFE, coupled on the surface GC=∂WFE, and projection of the surface coupling force onto the finite element using (b) the conventional and (c) CG-enriched approaches (adapted from [10]). 2.1 Governing equations 2.1.1 The microscopic particle model The Discrete Element Method models granular materials as assemblies of rigid particles that interact via binary inter-particle contact forces and torques. Simplified, rigid geometric objects (walls) can be added that also interact with the particles via binary contact forces and torques. External forces and torques can be applied to the particles as well. The macroscopic behavior of the particle assembly is then obtained by resolving the translational and rotational motions of all particles individually according to Newton’s second law. Consider particle athat has Nand Nwexisting contacts from the neighboring particles and walls. For given initial and boundary conditions, the time evolution of the velocity va, position xa, angular velocity wa, and orientation qaof a particle Pais given by dva dt=1 ma⇣N  b=1 fab + Nw  g=1 fw ag +fb a⌘ | {z } :=fa ,(1) dxa dt=va,(2) (3) dwa dt=I1 a⇣N  b=1 lab ⇥fab + Nw  g=1 lag ⇥fw ag⌘ | {z } :=ta ,(4) dqa dt=C(qa)wa.(5) where maand Iaare particle’s mass and inertia tensor, C(qa)is the transformation matrix that allows an efficient handling of its orientation qaas quaternion, fb aare the body forces, and faand taare the total forces and torques acting at the particle position xa. Eqs. (3) and (5) are closed using so-called contact laws that define inter-particle contact forces fab as a function of the overlap 1between a pair of interacting particles Paand Pb. Another contact law, similar to, or different from the inter-particle contact law, can be used to compute particle-wall interaction forces fw ag. The branch vector lab =xc ab xaconnects the particle center xawith the contact point xc ab, and the cross product lab ⇥fab contributes a torque that affects the particle rotation and angular velocity [12]. In its simplest form, a linear-elastic force-displacement law relates the contact force in the direction normal to the contact area and the inter-particle overlap un(resembling deformation) by a linear spring constant, namely, kfn abk=knun. Similar force-displacement laws are often taken in the tangential direction with a Coulomb yield criterion, kfs abkµkfn abk, in order to capture irrecoverable plastic deformation, using a friction coefficient µ. The above differential equations are solved numerically using the Velocity-Verlet algorithm. Using a higher-order accurate time integration scheme would not increase accuracy, because most contact laws used in the DEM have discontinuous derivatives. The open-source DEM code MercuryDPM [11] is used for the particle side of our new DEM-FEM coupled code. 2.1.2 The macroscopic continuum model We assume that the deformable body behaves as a three-dimensional elastic solid and describe its behavior using finite strain theory. A finite-strain framework is essential to correctly handle the impact of granular materials, in particular for deformable structures that have large displacements and/or geometrical nonlinearity. The vectors xand X, attached to the body, are used to define the material’s undeformed and deformed configurations. A strain measure, the Green strain tensor, is defined as e=1 2(FT·F1), where F =1+—u is the deformation gradient w.r.t. x,u=Xxthe displacement field, 1 the identity matrix, and (·)Tthe transpose. Let us consider a solid body subjected to a surface traction t on a subset Gtof the body’s boundary ∂WFE,a body force density b acting on the domain WFE, and a displacement boundary condition, prescribed by X(BC)on another subset of the boundary GX. Using the stress and strain measures sand e, the deformation is governed by 1Note that the incremental shear displacement computed from the relative velocity between two particles at the contact point is used to compute the tangential interaction force. Surface Coupling a smooth and continuous function w(x), x2⌦C,TW: You define w on X or x which110 monotonically increases from zero on @⌦DE \⌦Cto one on @⌦FE \⌦C, as shown in Fig. 1a.111 The virtual work done in the DEM model WDE is weighted by 1 w(x)suchthatthe112 coupling weights on the FEM and DEM sides sum to unity at any given location. The113 coupled governing equations also need to be adapted by the respective coupling weights, as114 demonstrated in Sec. 2.3.115 (a) (b) Figure 1: (a) 2D illustration of volume coupling between discrete particle and finite element domains ⌦FE and ⌦DE, in the overlapping volume ⌦C=⌦ FE \⌦DE =SNe e=1 ⌦C e, and (b) mapping of particle velocity v↵ onto FEM nodes. 2.2. Homogenization using coarse-graining116 Before applying constraints on the variables (e.g., displacement or velocity) at di↵erent117 scales, a computational homogenization is needed to bring them onto the same scale, as118 illustrated in Fig. 1b. The fact that FEM and DEM give solutions at the particle and119 continuum scales, respectively, makes their volume coupling non-trivial: the length scale120 (either microscopic or macroscopic) on which the variables should be constrained is not well121 defined. A novel approach we take here is CG-enriched homogenization, using smooth kernel122 functions whose length scale can be chosen independent of particle and finite element sizes.123 In the following, we will demonstrate how a computational homogenization is enriched with124 CG and in which limit this generalized formulation reduces to the conventional one.125 Let us now consider Npparticles that reside within a finite element in the coupling domain126 ⌦C e, as illustrated in Fig. 2. Using coarse-graining, we can define a homogenized velocity field127 6 Scale Coupling Timeline Tabletting Multi-resolved coupling Under-resolved simulations Change of plans Coupling particles to fluid Fully-resolved Multi-resolved Under-resolved ITrade-off between speed and detail of information. IFully-resolved: Lf1 8⇡dp,min IUnder-resolved: Lf3dp,max 8 Particle-Fluid Coupling 13 / 30
The Moomph framework – coupling discrete and continuum solvers Moomph coupling framework microdata macrodata hook hook Multiscale coupling Particle-fluid coupling Particle-solid coupling 5 124 5. ATHERMO-MECHANICALLY COUPLED MULTI-SCALE MODEL OF GRANULAR MEDIA Heat transfer model + sintering Figure 5.21: Multi-scale framework to simulate the laser sintering process using MercuryDPM and oomph-lib. solid finite element domain, there are a total of 1500 volume coupling elements. Each coupling element contains a single discrete particle. The dimensions of the domain have been chosen to be representative, based on the average size of PA12 particles, which is approximately 125 µm. For a visual representation of this configuration, see Figure 5.22. Figure 5.22: Coupled domain to simulate laser sintering. Discrete elements (particles) are modelled using MercuryDPM, finite elements are modelled using oomph-lib. The domains are coupled using volume elements, containing one discrete element per volume. The material parameters are listed in Table. 5.6 A temperature and pressure-dependent visco-elastic-plastic model [2] is utilized to 5.4. RESULTS AND DISCUSSION 5 117 Table 5.4: Parameters for simulation Parameter Units Value Length - ly[m] 25.0 Height - lz[m] 0.5 Width - lx[m] 0.5 Flexural rigidity - M/EI [N/m2] 1.0£10°38 Conductivity - Ø[W/mK] 1.0 Thermal expansion - Æ[1/±C] 0.0°7.1£10°3 Gravity - g[m/s2] 0.0 T|t=0=T0[±C] 0.0 T|z=ly=T18t>0[ ±C] 1.0 T|z=0=T28t>0[ ±C] 0.0 generalized Hooke’s law, and the Poisson effect is neglected. This corresponds to a steady heat transfer problem, for which is assumed that the thermal equilibrium is reached once the simulation starts. Two temperatures, T1=1.0 ±C and T2=0.0 ±C, are set to the top (z=lz) and bottom (z=0) surfaces of the domain, respectively. The beam is fully constrained at lx=0. Fig. 5.12 shows the cantilever beam and the thermal deflection. Figure 5.12: Thermal deflection of a cantilever beam. The beam of length lyis fully constrained at x=0. We assume no internal heat generation H=0, and no heat flux or heat transfer occurring across the boundaries of the volume V. Fig. 5.12 shows the cantilever beam deflecting. The displacement results for the nodes along the °yaxis, located at z=lz/2, with different thermal expansion coefficients Æ, are presented in Fig. 5.13. The cantilever beam deflects higher according to the increment of its thermal Multiphysics coupling 14 / 30
Particle-fluid coupling •Fully resolved: Accurate and well-defined, but computationally expensive. Valid for h < 1 8πdmin. •Under-resolved: Fast, but requires closure relations, modelling assumptions. Valid for h > 3dmax. 15 / 30
Fully resolved particle-fluid coupling 16 / 30
Under-resolved particle-fluid coupling We apply a two-phase flow model (Anderson & Jackson): •Each particle experiences a drag force, fd p=β(uf−up). •The fluid phase is governed by: ∂(ϵρ) ∂t +∇·(ϵρuf) = 0, ϵDuf Dt =−ϵ∇p+µ∇2uf+ϵg − fd. •Coupled via coarse-grained drag force and porosity (Rui & Heng): fd=X p fd pϕ(r −rp), ϵ =X p Vpϕ(r −rp) 17 / 30
Towards faster simulations Challenges •Mapping between particles and fluid mesh is expensive •Optimal coarse-graining mesh often differs from the fluid mesh Proposed Solution Use interpolation meshes that store mapping information and refine automatically Fluid solver Particle solver 18 / 30
Towards faster simulations Challenges •Mapping between particles and fluid mesh is expensive •Optimal coarse-graining mesh often differs from the fluid mesh Proposed Solution Use interpolation meshes that store mapping information and refine automatically Fluid solver Particle solver Interpolated fluid velocity Coarse-grained porosity/drag 18 / 30
Example: Particles settling in a lid-driven cavity Mesh Field Fluid velocity Porosity . 19 / 30
Validation: Particles settling in a lid-driven cavity ∼10 times speedup 20 / 30
Code development by MercuryLab This feature was developed by MercuryLab, a spin-off supporting MercuryDPM. What MercuryLab offers •Access to MercuryCloud (at-cost for academics) •Consultancy & training •Code support and maintenance Learn more at mercurylab.org or meet us at Powtech. All MercuryLab-developed features are contributed back to the open-source codebase. 21 / 30
What is ON-DEM Cost Action: Open Network on DEM Simulations •It is a COST-action •European funding for research networks •Funds: Meetings, workshops, visits, training, etc.. •300+ members •Open to academics and industry •Aims •promote extensive use of open-source DEM codes •produce scientific advances enabling DEM simulation of more realistic problems 26 / 30
ON-DEM Partners More ? Autopas 27 / 30
ON-DEM Partners More ? Autopas 27 / 30
ON-DEM Partners More ? Autopas 27 / 30
Join the Action! Cost Action: Open Network on DEM Simulations 1Passing through time & space scales GPU, multiscale, ML, coupled 2Going closer to physics Validation & calibration 3Data processing and visualisation Large data, common file formats 4Normalisation and best practices Benchmarking, validation, guidelines 5Industrial networking, commercial utilisation Training, UI/UX, Cloud computing There is a shortage of specialists: •Computer scientists •Applied mathematicians •Numerical analysts 28 / 30
ON-DEM Acknowledgements This talk and my participation in this conference was funded by the COST action ON-DEM Check out ON-DEM.org or scan the QR code to join 29 / 30
Conclusions and Outlook Conclusions •Surface, multi-scale and particle-fluid coupling implemented in Moomph and (coming) with preCICE. •Applied to various applications (wall impact, dosing into elastic pouch, thermal expansion, sedimentation). •Coarse-graining improves coupling by smoothing the coupled fields (less oscillations, improved energy conservation). 30 / 30
Conclusions and Outlook Conclusions •Surface, multi-scale and particle-fluid coupling implemented in Moomph and (coming) with preCICE. •Applied to various applications (wall impact, dosing into elastic pouch, thermal expansion, sedimentation). •Coarse-graining improves coupling by smoothing the coupled fields (less oscillations, improved energy conservation). Present and Future developments •Compare and validate particle-fluid coupling techniques. •Validate with experiments for rotating drum case. •Add preCICE adaptor to oomph-lib. •Include thermal and chemical effects. •Extend multi-scale modelling to process scale, using model order reduction. 30 / 30
Conclusions and Outlook Present and Future developments •Compare and validate particle-fluid coupling techniques. •Validate with experiments for rotating drum case. •Add preCICE adaptor to oomph-lib. •Include thermal and chemical effects. •Extend multi-scale modelling to process scale, using model order reduction. Broader Impact (open-source ecosystem for coupled simulations, ON-DEM) •ON-DEM is looking for wider members. •Dedicated to open-source. •Please consider joining. •Follow-up network on coupled problems is being written. 30 / 30
Choosing the right coarse-graining scale (Weinhart et al., 2013, Phys. Fl.) 1 / 39
Coarse graining functions Essential properties of ϕ: 1Mass conservation: V=RR3ϕ(r)dr= 1, 2Positivity: ϕ(r)≥0∀r∈R3, 3Locality: ϕ(r) = 0 for |r|> c. Commonly used functions: 1Cutoff Gaussian 2Heaviside 3Lucy polynomial |r| φ3 G φ3 ˜ G φ3 ˜ H φ3 ˜ L 0 1 23 0 0.2 0.4 CG functions of width w= 1 6 / 39
Coarse graining functions Essential properties of ϕ: 1Mass conservation: V=RR3ϕ(r)dr= 1, 2Positivity: ϕ(r)≥0∀r∈R3, 3Locality: ϕ(r) = 0 for |r|> c. Commonly used functions: 1Cutoff Gaussian 2Heaviside 3Lucy polynomial |r| φ3 G φ3 ˜ G φ3 ˜ H φ3 ˜ L 0 1 23 0 0.2 0.4 CG functions of width w= 1 6 / 39
Coarse graining: mass balance We define velocity us.t. mass balance, ∂ρ ∂t +∇·(ρu)=0, is satisfied: u=j ρ,j= n X i=1 miviϕ(r−ri). 7 / 39
Coarse graining: momentum balance Similarly, we define stress σbased on momentum balance: ∂ρu ∂t +∇·(ρuu) = −∇·σ+ρg. σ=σk+σc, σk=Xn i=1 miv′ iv′ iϕ(r−ri), σc=Xi,j fijbij Z1 0 ϕ(r−ri+sbij)ds 0 3 Magnitude of stress |σ|2 8 / 39
Advantages CG +Produces continuum fields (not discrete values) +Satisfies local mass/momentum conservation for all w,N.1 +General: Accounts for boundary interactions3, granular mixtures4, temporal smoothing4, non-spherical particles5, etc +Applicable to both DPM and experimental data6 +Available open-source at MercuryDPM.org +Used for two-way coupling, e.g. with oomph-lib7 1Goldhirsch, Granular Matter (2010) 12(3), 239–252 2Weinhart et al, Physics of Fluids (2013) 25(7), 070605 3Weinhart et al, Granular Matter (2012) 14(2), 289-294 4Tunuguntla et al, Comput. Particle Mechanics (2016) 3(3), 349-365 5Weinhart et al, Powder Technology (2016) 293, 138-148 6Roy et al, EPJ E (2019) 42: 14 7Cheng et al, CMAME, 403, 115651 (2023) 9 / 39
Continuum modelling with oomph-lib •Finite element solver developed at University of Manchester. •Available open-source at oomph-lib.maths.man.ac.uk •Focus on fluid-structure interactions. •Supports mesh refinement. •Similar code base as MercuryDPM 100 Comput Mech (2008) 43:91–101 Fig. 6 Flow past a cylinder with an attached flag for the parameter values of test case “FSI2” in Ref. [10]. aTime-trace of the x2-coordinate of the flag’s tip. bSnapshot of the flow field (instantaneous streamlines and pressure contours) at t=6.04 Fig. 7 Flow past a cylinder with an attached flag for the parameter values of test case “FSI3” in Ref. [10]. aTime-trace of the x2-coordinate of the flag’s tip. bSnapshot of the flow field (instantaneous streamlines and pressure contours) at t=3.615 time of 75.3 s for the monolithic Newton solver at each timestep. Figure 7shows the corresponding results for the test case “FSI3” where fluid and wall densities are equal and Re = 200,St =0.05,Q=3.571 ×10−6and Λ2=1.7855 × 10−6.Comparedtotheconditionsintestcase“FSI2”,the system performs oscillations of much higher frequency and smaller amplitude. Our computation was performed with a fixed timestep of ∆t=0.005 and resulted in oscillations with a period of ≈0.19 and an amplitude of the vertical tip displacementof≈0.01±0.36.Theincreasein frequencyand Reynolds number leads to the development of thinner boundary and shear layers which require a finer spatial resolution, involving an average of 84,000 degrees of freedom. This had very little effect on the GMRES convergence rates. With P P P1 used as preconditioner an average of 45.8 iterations were required to obtain a converged solution of the linear systems,whilethesolutionofthe nonlinearsystemsbyNewton’s method now required an average of 104.4 s. Fig. 8 Unsteady flow through a collapsible tube for Q=10−7,St = 1. aTime-trace of the radial position of a control point located at 70% of the tube length, at the most strongly collapsed point in the cross section, for various Reynolds numbers. bSnapshot of the axial velocity profiles at t=181.26 for the Re =100 case 4.2.4 Self-excited oscillations of a 3D collapsible tube In the final example we evaluate the performance of the FSI preconditioner in preliminary 3D simulations of unsteady finite-Reynolds number flows in a collapsible tube, modelled as a circular cylindrical shell of radius L,length10L and wall thickness h/L=1/20, and a Poisson ratio of ν=0.49, mounted on two short rigid tubes of length L. We prescribe Poiseuille flow with average velocity Uat the far downstream end and impose parallel, axially tractionfree inflow upstream. Figure 8shows the time-trace, rctrl(t), of the radial position of a control point on the tube wall, located towards the downstream end of the tube at 70% of its length, for Q=10−7,St =1andΛ=0atvarious Reynolds numbers. In all cases the simulation was started from the steady solution with rctrl =0.55, which required different external pressures for each different Reynolds number. At t=0wereducedpext to a value corresponding to an equilibrium configuration in which rctrl =0.6 and followed the system’s evolution. The time-trace shows that the system oscillates about the new equilibrium position. At small Reynolds numbers the oscillations decay but for sufficiently large Reynolds number they grow in amplitude, as observed in physical experiments; see e.g. [8]. Theresults presentedhere arepreliminary in thesense that the assessment of their mesh and timestep independence has notyetbeencompleted.However, theplot of thevelocityfield in Fig. 8suggests that even for a relatively coarse discretisation with 57,486 degrees of freedom, the velocity field is well resolved. The Jacobian matrix arising from the monolithic discretisation of 3D problems is significantly more dense 123 10 / 39
The Moomph framework – coupling discrete and continuum solvers Moomph coupling framework microdata macrodata hook hook Multiscale coupling Particle-fluid coupling Particle-solid coupling 5 124 5. ATHERMO-MECHANICALLY COUPLED MULTI-SCALE MODEL OF GRANULAR MEDIA Heat transfer model + sintering Figure 5.21: Multi-scale framework to simulate the laser sintering process using MercuryDPM and oomph-lib. solid finite element domain, there are a total of 1500 volume coupling elements. Each coupling element contains a single discrete particle. The dimensions of the domain have been chosen to be representative, based on the average size of PA12 particles, which is approximately 125 µm. For a visual representation of this configuration, see Figure 5.22. Figure 5.22: Coupled domain to simulate laser sintering. Discrete elements (particles) are modelled using MercuryDPM, finite elements are modelled using oomph-lib. The domains are coupled using volume elements, containing one discrete element per volume. The material parameters are listed in Table. 5.6 A temperature and pressure-dependent visco-elastic-plastic model [2] is utilized to 5.4. RESULTS AND DISCUSSION 5 117 Table 5.4: Parameters for simulation Parameter Units Value Length - ly[m] 25.0 Height - lz[m] 0.5 Width - lx[m] 0.5 Flexural rigidity - M/EI [N/m2] 1.0£10°38 Conductivity - Ø[W/mK] 1.0 Thermal expansion - Æ[1/±C] 0.0°7.1£10°3 Gravity - g[m/s2] 0.0 T|t=0=T0[±C] 0.0 T|z=ly=T18t>0[ ±C] 1.0 T|z=0=T28t>0[ ±C] 0.0 generalized Hooke’s law, and the Poisson effect is neglected. This corresponds to a steady heat transfer problem, for which is assumed that the thermal equilibrium is reached once the simulation starts. Two temperatures, T1=1.0 ±C and T2=0.0 ±C, are set to the top (z=lz) and bottom (z=0) surfaces of the domain, respectively. The beam is fully constrained at lx=0. Fig. 5.12 shows the cantilever beam and the thermal deflection. Figure 5.12: Thermal deflection of a cantilever beam. The beam of length lyis fully constrained at x=0. We assume no internal heat generation H=0, and no heat flux or heat transfer occurring across the boundaries of the volume V. Fig. 5.12 shows the cantilever beam deflecting. The displacement results for the nodes along the °yaxis, located at z=lz/2, with different thermal expansion coefficients Æ, are presented in Fig. 5.13. The cantilever beam deflects higher according to the increment of its thermal Multiphysics coupling 11 / 39
Particle-solid coupling (H Cheng) Concurrent multi-scale modeling of particle-continuum interactions 3 (a) (b) (c) Fig. 1: (a) 2D illustration of the discrete particle and finite element domains WDE and WFE, coupled on the surface GC=∂WFE, and projection of the surface coupling force onto the finite element using (b) the conventional and (c) CG-enriched approaches (adapted from [10]). 2.1 Governing equations 2.1.1 The microscopic particle model The Discrete Element Method models granular materials as assemblies of rigid particles that interact via binary inter-particle contact forces and torques. Simplified, rigid geometric objects (walls) can be added that also interact with the particles via binary contact forces and torques. External forces and torques can be applied to the particles as well. The macroscopic behavior of the particle assembly is then obtained by resolving the translational and rotational motions of all particles individually according to Newton’s second law. Consider particle athat has Nand Nwexisting contacts from the neighboring particles and walls. For given initial and boundary conditions, the time evolution of the velocity va, position xa, angular velocity wa, and orientation qaof a particle Pais given by dva dt=1 ma⇣N  b=1 fab + Nw  g=1 fw ag +fb a⌘ | {z } :=fa ,(1) dxa dt=va,(2) (3) dwa dt=I1 a⇣N  b=1 lab ⇥fab + Nw  g=1 lag ⇥fw ag⌘ | {z } :=ta ,(4) dqa dt=C(qa)wa.(5) where maand Iaare particle’s mass and inertia tensor, C(qa)is the transformation matrix that allows an efficient handling of its orientation qaas quaternion, fb aare the body forces, and faand taare the total forces and torques acting at the particle position xa. Eqs. (3) and (5) are closed using so-called contact laws that define inter-particle contact forces fab as a function of the overlap 1between a pair of interacting particles Paand Pb. Another contact law, similar to, or different from the inter-particle contact law, can be used to compute particle-wall interaction forces fw ag. The branch vector lab =xc ab xaconnects the particle center xawith the contact point xc ab, and the cross product lab ⇥fab contributes a torque that affects the particle rotation and angular velocity [12]. In its simplest form, a linear-elastic force-displacement law relates the contact force in the direction normal to the contact area and the inter-particle overlap un(resembling deformation) by a linear spring constant, namely, kfn abk=knun. Similar force-displacement laws are often taken in the tangential direction with a Coulomb yield criterion, kfs abkµkfn abk, in order to capture irrecoverable plastic deformation, using a friction coefficient µ. The above differential equations are solved numerically using the Velocity-Verlet algorithm. Using a higher-order accurate time integration scheme would not increase accuracy, because most contact laws used in the DEM have discontinuous derivatives. The open-source DEM code MercuryDPM [11] is used for the particle side of our new DEM-FEM coupled code. 2.1.2 The macroscopic continuum model We assume that the deformable body behaves as a three-dimensional elastic solid and describe its behavior using finite strain theory. A finite-strain framework is essential to correctly handle the impact of granular materials, in particular for deformable structures that have large displacements and/or geometrical nonlinearity. The vectors xand X, attached to the body, are used to define the material’s undeformed and deformed configurations. A strain measure, the Green strain tensor, is defined as e=1 2(FT·F1), where F =1+—u is the deformation gradient w.r.t. x,u=Xxthe displacement field, 1 the identity matrix, and (·)Tthe transpose. Let us consider a solid body subjected to a surface traction t on a subset Gtof the body’s boundary ∂WFE,a body force density b acting on the domain WFE, and a displacement boundary condition, prescribed by X(BC)on another subset of the boundary GX. Using the stress and strain measures sand e, the deformation is governed by 1Note that the incremental shear displacement computed from the relative velocity between two particles at the contact point is used to compute the tangential interaction force. (a) Simulates particle interactions with an elastic solid. (b) Requires coupling via surface traction field, t=P−fαγδ(x−xc αγ) (c) For upscaled particles, smooth surface traction: tCG =P−fαγϕ(x−xc αγ). Thus, solid responds to bulk rather than individual particles. Cheng et al, Computer Methods in Applied Mechanics and Engineering, 403, 115651 (2023) 12 / 39
Particle-solid coupling (H Cheng) Concurrent multi-scale modeling of particle-continuum interactions 3 (a) (b) (c) Fig. 1: (a) 2D illustration of the discrete particle and finite element domains WDE and WFE, coupled on the surface GC=∂WFE, and projection of the surface coupling force onto the finite element using (b) the conventional and (c) CG-enriched approaches (adapted from [10]). 2.1 Governing equations 2.1.1 The microscopic particle model The Discrete Element Method models granular materials as assemblies of rigid particles that interact via binary inter-particle contact forces and torques. Simplified, rigid geometric objects (walls) can be added that also interact with the particles via binary contact forces and torques. External forces and torques can be applied to the particles as well. The macroscopic behavior of the particle assembly is then obtained by resolving the translational and rotational motions of all particles individually according to Newton’s second law. Consider particle athat has Nand Nwexisting contacts from the neighboring particles and walls. For given initial and boundary conditions, the time evolution of the velocity va, position xa, angular velocity wa, and orientation qaof a particle Pais given by dva dt=1 ma⇣N  b=1 fab + Nw  g=1 fw ag +fb a⌘ | {z } :=fa ,(1) dxa dt=va,(2) (3) dwa dt=I1 a⇣N  b=1 lab ⇥fab + Nw  g=1 lag ⇥fw ag⌘ | {z } :=ta ,(4) dqa dt=C(qa)wa.(5) where maand Iaare particle’s mass and inertia tensor, C(qa)is the transformation matrix that allows an efficient handling of its orientation qaas quaternion, fb aare the body forces, and faand taare the total forces and torques acting at the particle position xa. Eqs. (3) and (5) are closed using so-called contact laws that define inter-particle contact forces fab as a function of the overlap 1between a pair of interacting particles Paand Pb. Another contact law, similar to, or different from the inter-particle contact law, can be used to compute particle-wall interaction forces fw ag. The branch vector lab =xc ab xaconnects the particle center xawith the contact point xc ab, and the cross product lab ⇥fab contributes a torque that affects the particle rotation and angular velocity [12]. In its simplest form, a linear-elastic force-displacement law relates the contact force in the direction normal to the contact area and the inter-particle overlap un(resembling deformation) by a linear spring constant, namely, kfn abk=knun. Similar force-displacement laws are often taken in the tangential direction with a Coulomb yield criterion, kfs abkµkfn abk, in order to capture irrecoverable plastic deformation, using a friction coefficient µ. The above differential equations are solved numerically using the Velocity-Verlet algorithm. Using a higher-order accurate time integration scheme would not increase accuracy, because most contact laws used in the DEM have discontinuous derivatives. The open-source DEM code MercuryDPM [11] is used for the particle side of our new DEM-FEM coupled code. 2.1.2 The macroscopic continuum model We assume that the deformable body behaves as a three-dimensional elastic solid and describe its behavior using finite strain theory. A finite-strain framework is essential to correctly handle the impact of granular materials, in particular for deformable structures that have large displacements and/or geometrical nonlinearity. The vectors xand X, attached to the body, are used to define the material’s undeformed and deformed configurations. A strain measure, the Green strain tensor, is defined as e=1 2(FT·F1), where F =1+—u is the deformation gradient w.r.t. x,u=Xxthe displacement field, 1 the identity matrix, and (·)Tthe transpose. Let us consider a solid body subjected to a surface traction t on a subset Gtof the body’s boundary ∂WFE,a body force density b acting on the domain WFE, and a displacement boundary condition, prescribed by X(BC)on another subset of the boundary GX. Using the stress and strain measures sand e, the deformation is governed by 1Note that the incremental shear displacement computed from the relative velocity between two particles at the contact point is used to compute the tangential interaction force. Concurrent multi-scale modeling of particle-continuum interactions 3 (a) (b) (c) Fig. 1: (a) 2D illustration of the discrete particle and finite element domains WDE and WFE, coupled on the surface GC=∂WFE, and projection of the surface coupling force onto the finite element using (b) the conventional and (c) CG-enriched approaches (adapted from [10]). 2.1 Governing equations 2.1.1 The microscopic particle model The Discrete Element Method models granular materials as assemblies of rigid particles that interact via binary inter-particle contact forces and torques. Simplified, rigid geometric objects (walls) can be added that also interact with the particles via binary contact forces and torques. External forces and torques can be applied to the particles as well. The macroscopic behavior of the particle assembly is then obtained by resolving the translational and rotational motions of all particles individually according to Newton’s second law. Consider particle athat has Nand Nwexisting contacts from the neighboring particles and walls. For given initial and boundary conditions, the time evolution of the velocity va, position xa, angular velocity wa, and orientation qaof a particle Pais given by dva dt=1 ma⇣N  b=1 fab + Nw  g=1 fw ag +fb a⌘ | {z } :=fa ,(1) dxa dt=va,(2) (3) dwa dt=I1 a⇣N  b=1 lab ⇥fab + Nw  g=1 lag ⇥fw ag⌘ | {z } :=ta ,(4) dqa dt=C(qa)wa.(5) where maand Iaare particle’s mass and inertia tensor, C(qa)is the transformation matrix that allows an efficient handling of its orientation qaas quaternion, fb aare the body forces, and faand taare the total forces and torques acting at the particle position xa. Eqs. (3) and (5) are closed using so-called contact laws that define inter-particle contact forces fab as a function of the overlap 1between a pair of interacting particles Paand Pb. Another contact law, similar to, or different from the inter-particle contact law, can be used to compute particle-wall interaction forces fw ag. The branch vector lab =xc ab xaconnects the particle center xawith the contact point xc ab, and the cross product lab ⇥fab contributes a torque that affects the particle rotation and angular velocity [12]. In its simplest form, a linear-elastic force-displacement law relates the contact force in the direction normal to the contact area and the inter-particle overlap un(resembling deformation) by a linear spring constant, namely, kfn abk=knun. Similar force-displacement laws are often taken in the tangential direction with a Coulomb yield criterion, kfs abkµkfn abk, in order to capture irrecoverable plastic deformation, using a friction coefficient µ. The above differential equations are solved numerically using the Velocity-Verlet algorithm. Using a higher-order accurate time integration scheme would not increase accuracy, because most contact laws used in the DEM have discontinuous derivatives. The open-source DEM code MercuryDPM [11] is used for the particle side of our new DEM-FEM coupled code. 2.1.2 The macroscopic continuum model We assume that the deformable body behaves as a three-dimensional elastic solid and describe its behavior using finite strain theory. A finite-strain framework is essential to correctly handle the impact of granular materials, in particular for deformable structures that have large displacements and/or geometrical nonlinearity. The vectors xand X, attached to the body, are used to define the material’s undeformed and deformed configurations. A strain measure, the Green strain tensor, is defined as e=1 2(FT·F1), where F =1+—u is the deformation gradient w.r.t. x,u=Xxthe displacement field, 1 the identity matrix, and (·)Tthe transpose. Let us consider a solid body subjected to a surface traction t on a subset Gtof the body’s boundary ∂WFE,a body force density b acting on the domain WFE, and a displacement boundary condition, prescribed by X(BC)on another subset of the boundary GX. Using the stress and strain measures sand e, the deformation is governed by 1Note that the incremental shear displacement computed from the relative velocity between two particles at the contact point is used to compute the tangential interaction force. (a) Simulates particle interactions with an elastic solid. (b) Requires coupling via surface traction field, t=P−fαγδ(x−xc αγ) (c) For upscaled particles, smooth surface traction: tCG =P−fαγϕ(x−xc αγ). Thus, solid responds to bulk rather than individual particles. Cheng et al, Computer Methods in Applied Mechanics and Engineering, 403, 115651 (2023) 12 / 39
Particle-solid coupling: Example Granular jet impacting an elastic beam 13 / 39
Multiscale coupling: Derivation Principle of virtual work: •Derive a continuum model by balancing the work done when applying infinitesimally small displacements δX: ZΩhσ:δϵ−(b−ρ¨ X)i·δXdV+ZΓt t·δXdA= 0 •Similarly, we can write the discrete model as Xαh(mα¨xα−fα)−Xβfαβi·δxα= 0 20 / 39
Multiscale coupling: Derivation Now couple these two approaches ZΩh(1 −w)σ:δϵ−(1 −w)(b−ρ¨ X)−fci·δXdV+ZΓt (1 −w)t·δXdA= 0 Xαhwα(mα¨xα−fα) + fc α−Xβwαβfαβi·δxα= 0 with coupling weights wand coupling forces fc=ϵ(vFE −vDE) fc α=ϵΠα(vFE −vDE) (a) v0 (b) Figure 18: (a) A regular granular square lattice with one third of its volume overlapping with a FEM mesh, and (b) injection of a compressional P-wave from the left end of the granular system at time t= 0s. The color indicates the coupling weight won the continuum body, varying from approximately 1.0 on the leftmost boundary to 0.0 on the rightmost boundary of the coupling zone. Note, one discrete particle per finite element is shown here only for the illustrative purpose; more particles per element are used for the simulations in Secs. 6.2.4–6.2.4. The initial velocity of the particles is indicated by the arrows in Figure 18b. Table 6: Penalty parameter for di↵erent finite particle-element configurations for volume coupling. Test cases 1P 2P 3P 4P Num. of particles per element n312 33343 Particle diameter d(mm) 400 200 133.¯ 3 100 Penalty parameter ✏(Pa/mm2)1⇥1044⇥1049⇥10416 ⇥104 Penalty parameter ✏d2(Pa) 4 ⇥108 40 Cheng et al, Computer Methods in Applied Mechanics and Engineering, 403, 115651 (2023) 21 / 39
Multiscale coupling: Derivation a smooth and continuous function w(x), x2⌦C,TW: You define w on X or x which110 monotonically increases from zero on @⌦DE \⌦Cto one on @⌦FE \⌦C, as shown in Fig. 1a.111 The virtual work done in the DEM model WDE is weighted by 1 w(x)suchthatthe112 coupling weights on the FEM and DEM sides sum to unity at any given location. The113 coupled governing equations also need to be adapted by the respective coupling weights, as114 demonstrated in Sec. 2.3.115 (a) (b) Figure 1: (a) 2D illustration of volume coupling between discrete particle and finite element domains ⌦FE and ⌦DE, in the overlapping volume ⌦C=⌦ FE \⌦DE =SNe e=1 ⌦C e, and (b) mapping of particle velocity v↵ onto FEM nodes. 2.2. Homogenization using coarse-graining116 Before applying constraints on the variables (e.g., displacement or velocity) at di↵erent117 scales, a computational homogenization is needed to bring them onto the same scale, as118 illustrated in Fig. 1b. The fact that FEM and DEM give solutions at the particle and119 continuum scales, respectively, makes their volume coupling non-trivial: the length scale120 (either microscopic or macroscopic) on which the variables should be constrained is not well121 defined. A novel approach we take here is CG-enriched homogenization, using smooth kernel122 functions whose length scale can be chosen independent of particle and finite element sizes.123 In the following, we will demonstrate how a computational homogenization is enriched with124 CG and in which limit this generalized formulation reduces to the conventional one.125 Let us now consider Npparticles that reside within a finite element in the coupling domain126 ⌦C e, as illustrated in Fig. 2. Using coarse-graining, we can define a homogenized velocity field127 6 a smooth and continuous function w(x), x2⌦C,TW: You define w on X or x which110 monotonically increases from zero on @⌦DE \⌦Cto one on @⌦FE \⌦C, as shown in Fig. 1a.111 The virtual work done in the DEM model WDE is weighted by 1 w(x)suchthatthe112 coupling weights on the FEM and DEM sides sum to unity at any given location. The113 coupled governing equations also need to be adapted by the respective coupling weights, as114 demonstrated in Sec. 2.3.115 (a) (b) Figure 1: (a) 2D illustration of volume coupling between discrete particle and finite element domains ⌦FE and ⌦DE, in the overlapping volume ⌦C=⌦ FE \⌦DE =SNe e=1 ⌦C e, and (b) mapping of particle velocity v↵ onto FEM nodes. 2.2. Homogenization using coarse-graining116 Before applying constraints on the variables (e.g., displacement or velocity) at di↵erent117 scales, a computational homogenization is needed to bring them onto the same scale, as118 illustrated in Fig. 1b. The fact that FEM and DEM give solutions at the particle and119 continuum scales, respectively, makes their volume coupling non-trivial: the length scale120 (either microscopic or macroscopic) on which the variables should be constrained is not well121 defined. A novel approach we take here is CG-enriched homogenization, using smooth kernel122 functions whose length scale can be chosen independent of particle and finite element sizes.123 In the following, we will demonstrate how a computational homogenization is enriched with124 CG and in which limit this generalized formulation reduces to the conventional one.125 Let us now consider Npparticles that reside within a finite element in the coupling domain126 ⌦C e, as illustrated in Fig. 2. Using coarse-graining, we can define a homogenized velocity field127 6 (a) Coupling force, fc=ϵ(vFE −vDE), penalises velocity differences. (b) CG used to produce a smooth velocity field, vDE =1 ρPimiviϕ(r−ri). 22 / 39
Multiscale coupling: Derivation a smooth and continuous function w(x), x2⌦C,TW: You define w on X or x which110 monotonically increases from zero on @⌦DE \⌦Cto one on @⌦FE \⌦C, as shown in Fig. 1a.111 The virtual work done in the DEM model WDE is weighted by 1 w(x)suchthatthe112 coupling weights on the FEM and DEM sides sum to unity at any given location. The113 coupled governing equations also need to be adapted by the respective coupling weights, as114 demonstrated in Sec. 2.3.115 (a) (b) Figure 1: (a) 2D illustration of volume coupling between discrete particle and finite element domains ⌦FE and ⌦DE, in the overlapping volume ⌦C=⌦ FE \⌦DE =SNe e=1 ⌦C e, and (b) mapping of particle velocity v↵ onto FEM nodes. 2.2. Homogenization using coarse-graining116 Before applying constraints on the variables (e.g., displacement or velocity) at di↵erent117 scales, a computational homogenization is needed to bring them onto the same scale, as118 illustrated in Fig. 1b. The fact that FEM and DEM give solutions at the particle and119 continuum scales, respectively, makes their volume coupling non-trivial: the length scale120 (either microscopic or macroscopic) on which the variables should be constrained is not well121 defined. A novel approach we take here is CG-enriched homogenization, using smooth kernel122 functions whose length scale can be chosen independent of particle and finite element sizes.123 In the following, we will demonstrate how a computational homogenization is enriched with124 CG and in which limit this generalized formulation reduces to the conventional one.125 Let us now consider Npparticles that reside within a finite element in the coupling domain126 ⌦C e, as illustrated in Fig. 2. Using coarse-graining, we can define a homogenized velocity field127 6 a smooth and continuous function w(x), x2⌦C,TW: You define w on X or x which110 monotonically increases from zero on @⌦DE \⌦Cto one on @⌦FE \⌦C, as shown in Fig. 1a.111 The virtual work done in the DEM model WDE is weighted by 1 w(x)suchthatthe112 coupling weights on the FEM and DEM sides sum to unity at any given location. The113 coupled governing equations also need to be adapted by the respective coupling weights, as114 demonstrated in Sec. 2.3.115 (a) (b) Figure 1: (a) 2D illustration of volume coupling between discrete particle and finite element domains ⌦FE and ⌦DE, in the overlapping volume ⌦C=⌦ FE \⌦DE =SNe e=1 ⌦C e, and (b) mapping of particle velocity v↵ onto FEM nodes. 2.2. Homogenization using coarse-graining116 Before applying constraints on the variables (e.g., displacement or velocity) at di↵erent117 scales, a computational homogenization is needed to bring them onto the same scale, as118 illustrated in Fig. 1b. The fact that FEM and DEM give solutions at the particle and119 continuum scales, respectively, makes their volume coupling non-trivial: the length scale120 (either microscopic or macroscopic) on which the variables should be constrained is not well121 defined. A novel approach we take here is CG-enriched homogenization, using smooth kernel122 functions whose length scale can be chosen independent of particle and finite element sizes.123 In the following, we will demonstrate how a computational homogenization is enriched with124 CG and in which limit this generalized formulation reduces to the conventional one.125 Let us now consider Npparticles that reside within a finite element in the coupling domain126 ⌦C e, as illustrated in Fig. 2. Using coarse-graining, we can define a homogenized velocity field127 6 (a) Coupling force, fc=ϵ(vFE −vDE), penalises velocity differences. (b) CG used to produce a smooth velocity field, vDE =1 ρPimiviϕ(r−ri). 22 / 39
Multiscale coupling: Wave propagation in elastic beam (a) Total mass (b) Sub-model linear momentum (c) Total linear momentum (d) Total energy Figure 20: Evolution of (a) mass, (b-c) momentum, and (d) kinetic, (e) potential and (f) total energies in the coupled FEM-DEM model, using the conventional and CG-enriched formulations (with various coarsegraining widths). instability and a converged solution cannot be obtained. This is because when the particle691 velocity is homogenized over too many finite elements, the di↵erence between uDE and uFE 692 is excessively large, especially for those elements which do not overlap with the particles.693 As a result, the coupling terms that arise by penalizing uDE uFE lead to non-converging694 iterations in finding the solution to approach uDE =uFE. As shown in Fig. 21, for a sharp695 signal, the CG width should not be larger than twice the element size.696 6.2.3. From a sharp impulse to continuous waveforms697 Although the previously considered boundary condition for the DEM model (i.e., free698 end on the left in Fig. 18b) ensures exact mass, momentum, and energy conservation, it699 43 Coarse-graining reduces energy loss in overlap region. 23 / 39
In development: Thermo-mechanical (Juan Alvarez) Apply a similar coupling to thermomechanical problems (elasticity + heat), using the principles of virtual heat. The number of elements for the FEM analysis is 4 for xand z,and50 for ydirection. The constitutive relation of the beam is represented by the generalized Hooke’s law, and the Poisson e↵ect is not taken into account. Two di↵erent temperatures T1=1.0CandT2Caresettothetopand bottom surface of the domain, respectively. Fig. 13 shows the deflection of a cantilever beam under di↵erent ↵. Figure 13: Thermal deflection results. Fig. 13 shows di↵erent deflections when increasing the thermal inertia coefficient (↵). The comparison with the analytical solution is presented in Fig. 14 where the deflection of the cantilever beam for di↵erent thermal inertia values is shown. The agreement between the analytical solution and simulations is consistent throughout the range of ↵=0.07.1⇥1031/C. This indicates that the continuum analysis using isotropic material properties may be reliable for predicting thermal deflection. 42 The number of elements for the FEM analysis is 4 for xand z,and50 for ydirection. The constitutive relation of the beam is represented by the generalized Hooke’s law, and the Poisson e↵ect is not taken into account. Two di↵erent temperatures T1=1.0CandT2Caresettothetopand bottom surface of the domain, respectively. Fig. 13 shows the deflection of a cantilever beam under di↵erent ↵. Figure 13: Thermal deflection results. Fig. 13 shows di↵erent deflections when increasing the thermal inertia coefficient (↵). The comparison with the analytical solution is presented in Fig. 14 where the deflection of the cantilever beam for di↵erent thermal inertia values is shown. The agreement between the analytical solution and simulations is consistent throughout the range of ↵=0.07.1⇥1031/C. This indicates that the continuum analysis using isotropic material properties may be reliable for predicting thermal deflection. 42 Bending of an elastic beam due to heat. 24 / 39
Particle-fluid coupling •Fully resolved: Accurate and well-defined, but computationally expensive. Valid for h < 1 8πdmin. •Under-resolved: Fast, but requires closure relations, modelling assumptions. Valid for h > 3dmax. 25 / 39
Fully resolved particle-fluid coupling 26 / 39
Under-resolved particle-fluid coupling We apply a two-phase flow model (Anderson & Jackson): •Each particle experiences a drag force, fd p=β(uf−up). •The fluid phase is governed by: ∂(ϵρ) ∂t +∇·(ϵρuf) = 0, ϵDuf Dt =−ϵ∇p+µ∇2uf+ϵg−fd. •Coupled via coarse-grained drag force and porosity (Rui & Heng): fd=X p fd pϕ(r−rp), ϵ =X p Vpϕ(r−rp) 27 / 39
Underresolved particle-fluid coupling Particles settling in a lid-driven cavity 28 / 39
Open-source coupling with 34 / 39
Coupling with preCICE (with David Schneider) •Allows coupling between multiple codes and models via continuum fields. •Takes care of coupling details (mapping, time-step interpolation, etc) •Here, we will use MercuryDPM and openFoam (multiphaseEuler): Adapter ϵ(r), fd(r) uf(r) Continuum model 35 / 39
Coupling with preCICE (with David Schneider) •Allows coupling between multiple codes and models via continuum fields. •Takes care of coupling details (mapping, time-step interpolation, etc) •Here, we will use MercuryDPM and openFoam (multiphaseEuler): Adapter ϵ(r), fd(r) uf(r) Continuum model 35 / 39
Validation: Oscillating flows Rayleigh-Kothe flows are standard CFD/FSI test cases to test solver accuracy. Here we use it to compare mapping functions (nearest-neighbour vs rbf): •One-directional coupling (fluid →particles). •No inertia (particles follow streamlines). •Oscillating flow (particles return to original location). 36 / 39
Validation: Oscillating flows •Observation: rbf produces a smoother velocity field, reversible. 37 / 39
Validation: Sedimenting particles •Bidirectional coupling •Result: steady-state predicted well •But: dynamics sensitive to mapping -0.035 -0.03 -0.025 -0.02 -0.015 -0.01 -0.005 0 0 0.05 0.1 0.15 0.2 0.25 w = 5x10-3 w = 4x10-3 empirical correlation Velocity Time pressure velocity magnitude 38 / 39
Validation: Sedimenting particles •Bidirectional coupling •Result: steady-state predicted well •But: dynamics sensitive to mapping -0.035 -0.03 -0.025 -0.02 -0.015 -0.01 -0.005 0 0 0.05 0.1 0.15 0.2 0.25 w = 5x10-3 w = 4x10-3 empirical correlation Velocity Time pressure velocity magnitude 38 / 39
Coupling with preCICE: Next steps •Compare with Moomph coupling. •Validate with experiments for rotating drum case. •Use coarse-graining to extract coupled fields. 39 / 39
Goldhirsch, I. 2010 Stress, stress asymmetry and couple stress: from discrete particles to continuous fields. Granular Matter 12 (3), 239–252. Luding, S. 2008 Introduction to discrete element methods DEM: Basics of contact force models and how to perform the micro-macro transition to continuum theory. Euro. J. of Enviro. Civ. Eng. 12 (7-8), 785–826. Roy, S., Scheper, B. J., Polman, H., Thornton, A. R., Tunuguntla, D. R., Luding, S. & Weinhart, T. 2019 Surface floow profles for dry and wet granular materials by particle . Thornton, A. R. & Weinhart, T. et al. 2009-2019 Mercurydpm. http://MercuryDPM.org/. Tunuguntla, D. R., Weinhart, T. & Thornton, A. R. 2017 Comparing and contrasting size-based particle segregation models. Computational Particle Mechanics 4(4), 387–405. Weinhart, T., Hartkamp, R., Thornton, A. R. & Luding, S. 2013 Coarse-grained local and objective continuum description of 3d granular flows down an inclined surface. Phys. Fluids 25 (070605). 39 / 39
Weinhart, T., Thornton, A. R., Luding, S. & Bokhove, O. 2012 From discrete particles to continuum fields near a boundary. Granular Matter 14, 289–294. 39 / 39