Numerical Modelling of Granular Flows B. Maury Laboratoire de Math - PowerPoint PPT Presentation

Numerical Modelling of Granular Flows B. Maury Laboratoire de Math´ ematiques, Universit´ e Paris-Sud, Orsay Model (hard spheres, non elastic collisions) Strategies A stable method to handle non elastic collisions Numerical tests Theoretical framework and analysis of the scheme Gluey particle model (with A. Lefebvre)

R d (with d = 3, 2, or 1), radii ( r i ) 1 ≤ i ≤ N , mass 1 . R dN } N rigid spheres in Q = { q = ( q 1 , q 2 , · · · , q N ) ∈ Feasible set : Q 0 = { q ∈ Q, D ij ( q ) ≥ 0 ∀ i, j , 1 ≤ i < j ≤ N } where D ij ( q ) = | q i − q j | − ( r i + r j ) (Rk. : Obstacle treated in the same way: D ik ( q ) ≥ 0) N.B.: D ij may be negative, and it is smooth in a neighbourhood of Q 0 (for finite size particles).

e ij r j q j D ij r i q i R dN : gradient of the distance between spheres i and − e ij . G ij = ∇ D ij ∈ T Q = j : G ij = ( · · · , 0 , − e ij , 0 , · · · , 0 , e ij , 0 , · · · ) , e ij = q j − q i | q j − q i | ,

C q : set of feasible directions at q ∈ Q 0 , C q = { v ∈ T Q , G ij · v ≥ 0 as soon as D ij ( q ) = 0 } , Outward normal cone to Q 0 at q (polar cone of C q ): C ◦ = q = { h ∈ T Q , h · v ≤ 0 ∀ v ∈ C q } N q � = {− µ ij G ij ( q ) , µ ij ≥ 0 , D ij ( q ) µ ij = 0 } (Farkas Lemma) i<j C q U I d = P C q + P N q (Moreau) N q

→ q ( t ) ∈ Q 0 , u = ˙ Find t �− q ,  d u � = f ( q , t ) + λ ij G ij    dt   i<j   λ ij ∈ M + ( I ) , supp( λ ij ) ⊂ { t , D ij ( q ( t )) = 0 } ,       u + P C q u − ,  = C q = { v ∈ T Q , G ij · v ≥ 0 as soon as D ij ( q ) = 0 } , Extended collision model (with e ∈ [0 , 1]: restitution coefficient): u + = u − − (1 + e ) P N q u − , Differential inclusion formulation : d 2 q dt 2 + N q ∋ f

1D example q 1 ] τ, + ∞ [ − g λ = gτδ τ + g q t

1 D problems : Q 0 defined as a connected component of the set of all configurations q with no overlapping. R N , q i +1 − q i ≥ r i +1 + r i } In this case, Q 0 = { q = ( q 1 , . . . , q N ) ∈ is closed and convex, and N q is with the subdifferential of the indicatrix of Q 0 : � � 0 if q ∈ Q 0 � N q = ∂I Q 0 ( q ) with I Q 0 ( q ) = � � + ∞ if q / ∈ Q 0 � and hence q �− → N q is a maximal monotone operator. In general, Q 0 is prox-regular (see Federer [10], Edmont-Thibault [8]).

Theoretical analysis Schatzmann [22], Ballard [1], Moreau, Buttazzo [5]. Only analiticity ensures uniqueness. Numerical simulation Molecular Dynamics: slight overlapping allowed, short-range repulsive force with local damping (Glowinski [11], Luding [13], Richefeu). Contact Dynamics: (1) Prediction of the violated constraints, then succession of single contact problems, with relaxation (Moreau and coworkers in Montpellier [21]) (2) Linearization of the constraints, global handling of contacts (Stewart [23], Maury [17]).

� t n +1 f n +1 ( q ) = f ( q , t ). Numerical Scheme h t n 1. Initialization ( q 0 h , u 0 h ) = ( q 0 , u 0 ) . 2. Compute u n +1 as the solution to the constrained minimization h problem 1 � � 2 � � u − u n h − h f n +1 ( q n min h ) h 2 u ∈ C h ( q n h ) with C h ( q n h ) = { u ∈ T Q , D ij ( q n h ) + h G ij ( q n h ) · u ≥ 0 } . � �� � ≈ D ij ( q n h + h u ) 3. Update the positions q n +1 h + h u n +1 = q n . h h

Interpretation of the scheme u n +1 h ) ( u n h + h f n +1 ( q n = P C h ( q n h )) h h equivalent to say that u n h + h f n +1 ( q n h ) − u n +1 h ) ( u n +1 ∈ ∂I C h ( q n ) . h h h where � � 0 if v ∈ K � ∂ϕ ( x ) = { v, ϕ ( x )+( v, h ) ≤ ϕ ( x + h ) ∀ h } , I K ( v ) = � � + ∞ if v / ∈ K �

D 13 < 0 q 1 D 34 < 0 D 12 < 0 C q 1 q 2 C q 2 As a consequence, the scheme can be written u n +1 − u n � � u n +1 ∋ f n +1 ( q n h h + ∂I C h ( q n h ) h ) h h h

For any q ∈ Q 0 , N q is the convex closure of half lines − R + G ij for u n +1 − u n � � u n +1 ∋ f n +1 h h ( q n + ∂I C h ( q n h ) ( ⋆ ) h ) h h h indices verifying D ij ( q ) = 0, and ∂I C h ( q ) ( u ) can be written as the same sum for indices i and j such that D ij ( q ) + h G ij · u = 0 . Left-hand side: Taylor expansion of D ij ( q + h u ) h ) ( u n +1 − → the set ∂I C h ( q n ) can be seen as a prediction of N q n +1 . h h To sum up : ( ⋆ ) is a semi-implicit time discretization of inclusion d u dt + N q ∋ f ( q ) .

Interpretation in terms of positions forbidden forbidden forbidden

Projection step (independent from the scheme itself) with B ∈ M r, 2 N ( R ) ( r = N ( N − 1) / 2 = number of constraints). 1 2 | u − U | 2 , K = { v , B v ≤ D } Find u = arg min K

R 2 N × R r Dual formulation: Find ( u , λ ) ∈ + saddle point of R r L ( v , µ ) = 1 2 | v − U | 2 + ( B v − D , µ ) , L ( u , µ ) ≤ L ( u , λ ) ≤ L ( v , λ ) , ∀ v ∈ T Q , µ ∈ + . Equivalently:  u + B ⋆ λ = U      B u ≤ D      ( B u − D , λ ) = 0 . with U predicted velocity, u actual velocity. Uzawa algorithm: � � �� λ k +1 = Π + λ k + ρ B ( U − B ⋆ λ k ) − D

Back to our problem: U = u n h + h f n +1 ( q n h ) , B v = ( − h G ij · v ) i<j , D = (D ij ( q n h )) i<j , h so that u n +1 and λ n +1 = ( λ n +1 ) 1 ≤ i<j ≤ r are related by h h ij u n +1 − u n � = f n +1 λ n +1 h h ( q n G ij ( q n h ) + h ) . h ij h i<j d u � dt = f + λ ij G ij i<j

Behaviour of uzawa algorithm

Cost reduction Number of constraints r = N ( N − 1) / 2 can be reduced to O ( N ). Bucket sorting:

NUMERICAL TESTS Behaviour of the scheme in the case of huge time steps Confrontation to a case of non-uniqueness Many-body problem with large time steps Stochastic forcing Various animations

1D problem, spheres in a row, non elastic chock Discrete version of pressureless gas models See Brenier [4] (sticky particles) for punctual particles ∂ t ρ + ∂ x ( ρu ) = 0 , ∂ t ( ρu ) + ∂ x ( ρu 2 ) = 0 . or Berthelin [2] (sticky blocks) for finite size particule ∂ t ρ + ∂ x ( ρu ) = 0 , ∂ t ( ρu ) + ∂ x ( ρu 2 ) + ∂ x p = 0 . ρ ≤ 1 (1 − ρ ) p = 0 .

1.0 1.0 0.9 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0.0 0.0 0.0 0.6 1.2 1.8 2.4 3.0 0.0 0.6 1.2 1.8 2.4 3.0 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0

Non uniqueness (counter example in the spirit of Schatzman [22], Ballard [1]) u (0) = 0 , � t q ( t ) = u ( s ) ds ∀ t ∈ I, 0 u ( t ) ˙ = f ( t ) + λ, R − as soon as q = 0. supp( λ ) ⊂ { t , q ( t ) = 0 } , u − − P N q u − u + = ∀ t ∈ I, where N q is { 0 } whenever q > 0, and Z , k ≥ − 4 . Force field � � �  1 1 1 �  1 for t ∈ 2 k +1 , 2 k +1 + �   2 k +2 �  � f ( t ) = k ∈ � � � 1 2 k +2 , 1 1 �   − α for t ∈ 2 k +1 + �   2 k �

1 1 0.5 0.5 0 0 -0.5 -0.5 -1 -1 -1.5 -1.5 -2 -2 1/4 1/4 1/2 1/2 1 1 2 2 4 4

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1/4 1/2 1 2 4 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1/4 1/2 1 2 4

Large time steps in a many-body situation (every 3 time steps shots)

Stochastic forcing � t q ( t ) = q 0 + u ( s ) ds ∈ Q 0 ∀ t ∈ I 0 � d u = f ( q , t ) dt + dλ ij G ij ( q ( t )) + σd w i<j supp( dλ ij ) ⊂ { t, D ij ( q ( t )) = 0 } , u + ( t ) = P C q u − ( t ) ∀ t ∈ I, Numerically: √ � h σ w n +1 � 1 u n +1 h − h f n +1 � u − u n ( q n � � = arg min h ) − � h h 2 u ∈ C h ( q n h )

1D problem Unconstrained problem : primitive of the Brownian motion � t q = 0 w ( s ) ds . { t , q ( t ) = 0 } has a.s. a single cluster point at 0, Constrained problem : du = σdw + dλ dq = u dt u + P C q u − = Almost surely: points of Z = { t , q ( t ) = 0 } are left-hand cluster points of Z itself, but for a countable infinity (take off instants). N.B. : u is no longer in BV. Analysis : see Bertoin [3]

0.12 0.0032 0.0028 0.10 0.0024 0.08 0.0020 0.06 0.0016 0.0012 0.04 0.0008 0.02 0.0004 0.00 0.0000 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.00 0.01 0.02 0.03 0.04 0.05 9e-05 1e-06 9e-07 8e-05 8e-07 7e-05 7e-07 6e-05 6e-07 5e-05 5e-07 4e-05 4e-07 3e-05 3e-07 2e-05 2e-07 1e-05 1e-07 0e+00 0e+00 0.000 0.001 0.002 0.003 0.004 0.005 0.0e+00 4.0e-05 8.0e-05 1.2e-04 1.6e-04 2.0e-04 2.4e-04 2.8e-04 3.2e-04 3.6e-04 4.0e-04

Capillary forces (with common radius r ): F ji = − F ij = πγr (exp ( − AD ij + B ) + C ) e ij if D ij < D rupt . → 41-42

Further experiments Many-body simulations → 50-52-54 Macroscopic bodies Additional force : F = −∇ q V , � q i +1 − q i V = k ( | q i +1 − q i | − ℓ ) 2 + k a | q i +1 − q i | · q i − 1 − q i � | q i − 1 − q i | . 2 2 i − 1 i + 1 i → 62-64

First-order evolution equation : a crowd motion model (with Juliette Venel [16]) .

Safe zone ω Evolution equation : d q d q dt = P C q U ( q ) or dt + N q ∋ U . → 70-72-73