Non-homogeneous incompressible Bingham flows with variable yield - - PowerPoint PPT Presentation

Non-homogeneous incompressible Bingham flows with variable yield stress and application to volcanology.

Jordane Mathé with Laurent Chupin (maths) and Karim Kelfoun (LMV)

Laboratoire de Maths & Laboratoire Magmas et Volcans

june 2015

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 1 / 18

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 2 / 18

Laboratory experiment

Zoom into the granular column →

• O. Roche, LMV.

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 3 / 18

Laboratory experiment

Zoom into the granular column → Fluidisation : injection of gas through a pore plate.

• O. Roche, LMV.

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 3 / 18

Figure: Non-fluidised → short runout distance Figure: Fluidised → long runout distance

1

Model

2

Numerical simulation

3

Perspectives

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 5 / 18

1

Model

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 6 / 18

Two phases for one fluid

Consider one mixed fluid with variable density.

Two phases for one fluid

Consider one mixed fluid with variable density. It satisfies the non-homogeneous incompressible Navier-Stokes equations:

      

div(v) = 0 ∂tρ + div(ρv) = 0 ∂t(ρv) + div(ρv ⊗ v) + ∇p = ρg + div(S)

Two phases for one fluid

Consider one mixed fluid with variable density. It satisfies the non-homogeneous incompressible Navier-Stokes equations:

      

div(v) = 0 ∂tρ + div(ρv) = 0 ∂t(ρv) + div(ρv ⊗ v) + ∇p = ρg + div(S) Unknown: v: velocity, p: total pressure, ρ: density.

Two phases for one fluid

Consider one mixed fluid with variable density. It satisfies the non-homogeneous incompressible Navier-Stokes equations:

      

div(v) = 0 ∂tρ + div(ρv) = 0 ∂t(ρv) + div(ρv ⊗ v) + ∇p = ρg + div(S) Constant g: gravity.

Two phases for one fluid

Consider one mixed fluid with variable density. It satisfies the non-homogeneous incompressible Navier-Stokes equations:

      

div(v) = 0 ∂tρ + div(ρv) = 0 ∂t(ρv) + div(ρv ⊗ v) + ∇p = ρg + div(S) Constant g: gravity.

Rheology

Let precise S.

Rheology

In our case S = µDv + where Dv = ˙ γ is the strain rate tensor, µ is the effective viscosity

Rheology

In our case S = µDv + Σ , where Dv = ˙ γ is the strain rate tensor, µ is the effective viscosity Σ =q Dv |Dv| ⇒ Bingham fluid with yield stress q

Rheology

In our case S = µDv + Σ , where Dv = ˙ γ is the strain rate tensor, µ is the effective viscosity Σ =q Dv |Dv| ⇒ Bingham fluid with yield stress q

Idea

Let vary the yield stress q as a function of the interstitial gas pressure.

Variation of the yield stress

Definition of the yield stress

q =

• atmospheric pressure:

Coulomb friction high pressure: fluid

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 9 / 18

Variation of the yield stress

Definition of the yield stress

q =

• atmospheric pressure:

Coulomb friction high pressure: fluid q =

• tan(δ)(ρgh − gas pressure)

if low gas pressure if high gas pressure where δ is the internal friction angle.

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 9 / 18

Variation of the yield stress

Definition of the yield stress

q =

• atmospheric pressure:

Coulomb friction high pressure: fluid q =

• tan(δ)(ρgh − gas pressure)

if low gas pressure if high gas pressure q = tan(δ)(ρgh − gas pressure)+ where δ is the internal friction angle.

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 9 / 18

Equations of the model

Noting pf the interstitial gas pressure, we obtain:

                

div(v) = 0 ∂tρ + div(ρv) = 0 ∂t(ρv) + div(ρv ⊗ v) − µ∆v + ∇p = tan(δ) div

 (ρgh−pf )+

• yield q

Dv |Dv|

  + ρg

Equations of the model

Noting pf the interstitial gas pressure, we obtain:

                

div(v) = 0 ∂tρ + div(ρv) = 0 ∂t(ρv) + div(ρv ⊗ v) − µ∆v + ∇p = tan(δ) div

 (ρgh−pf )+

• yield q

Dv |Dv|

  + ρg

∂tpf + v · ∇pf − κ∆pf = 0 where κ is the diffusion coefficient.

2

Numerical simulation

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 11 / 18

Dambreak: Mesh

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 12 / 18

Dambreak: how to treat the density?

∂tρ + v · ∇ρ = 0

numerical method:

RK3 TVD - WENO5 scheme.

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 12 / 18

Numerical scheme (without density)

n = 0 v0, Σ0, p0 and pf 0 given.

Numerical scheme (without density)

n = 0 v0, Σ0, p0 and pf 0 given. n 0 vn, Σn, pn and pf n being calculated.

Numerical scheme (without density)

n = 0 v0, Σ0, p0 and pf 0 given. n 0 vn, Σn, pn and pf n being calculated. We compute ( vn+1, Σn+1) solution of

          

3 v n+1 − 4v n + v n−1 2δt + 2v n·∇v n − v n−1·∇v n−1 − ∆ v n+1 ρn+1 + ∇pn ρn+1 = div Σn+1 ρn+1 − ey, Σn+1 = Pqn(Σn+1 + r D v n+1 + ε(Σn − Σn+1)),

• v n+1
• ∂Ω

= 0.

Numerical scheme (without density)

n = 0 v0, Σ0, p0 and pf 0 given. n 0 vn, Σn, pn and pf n being calculated. We compute ( vn+1, Σn+1) solution of

          

3 v n+1 − 4v n + v n−1 2δt + 2v n·∇v n − v n−1·∇v n−1 − ∆ v n+1 ρn+1 + ∇pn ρn+1 = div Σn+1 ρn+1 − ey, Σn+1 = Pqn(Σn+1 + r D v n+1 + ε(Σn − Σn+1)),

• v n+1
• ∂Ω

= 0.

We compute (vn+1, pn+1) thanks to the incompressibility constrain

     v n+1 − v n+1 δt + 2 3ρn+1 ∇(pn+1 − pn) = 0, div v n+1 = 0, v n+1 · n

• ∂Ω = 0.

Numerical scheme (without density)

n = 0 v0, Σ0, p0 and pf 0 given. n 0 vn, Σn, pn and pf n being calculated. We compute ( vn+1, Σn+1) solution of

          

3 v n+1 − 4v n + v n−1 2δt + 2v n·∇v n − v n−1·∇v n−1 − ∆ v n+1 ρn+1 + ∇pn ρn+1 = div Σn+1 ρn+1 − ey, Σn+1 = Pqn(Σn+1 + r D v n+1 + ε(Σn − Σn+1)),

• v n+1
• ∂Ω

= 0.

We compute (vn+1, pn+1) thanks to the incompressibility constrain

     v n+1 − v n+1 δt + 2 3ρn+1 ∇(pn+1 − pn) = 0, div v n+1 = 0, v n+1 · n

• ∂Ω = 0.

We compute pf n+1 solution of

3pf n+1 − 4pf n + pf n−1 2δt + 2v n+1 · ∇pf

n − v n+1 · ∇pf n−1 − ∆pf n+1 = 0 .

Numerical scheme (without density)

n = 0 v0, Σ0, p0 and pf 0 given. n 0 vn, Σn, pn and pf n being calculated. We compute ( vn+1, Σn+1) solution of

          

3 v n+1 − 4v n + v n−1 2δt + 2v n·∇v n − v n−1·∇v n−1 − ∆ v n+1 ρn+1 + ∇pn ρn+1 = div Σn+1 ρn+1 − ey, Σn+1 = Pqn(Σn+1 + r D v n+1 + ε(Σn − Σn+1)),

• v n+1
• ∂Ω

= 0.

֒ → k = 0 Σn,0 = Σn, and (vn, pn, pn

f ) given.

֒ → k 0 Σn,k known, we first compute vn,k solution of a Laplace-type problem then we project the stress tensor to obtain Σn,k+1:

          

3 v n,k − 4v n + v n−1 2δt + 2v n·∇v n − v n−1·∇v n−1 − ∆ v n,k ρn+1 + ∇pn ρn+1 = div Σn,k ρn+1 − ey,

• v n,k
• ∂Ω

= 0, Σn,k+1 = Pqn(Σn,k + r D v n,k + ε(Σn − Σn,k)).

Numerical scheme (without density)

n = 0 v0, Σ0, p0 and pf 0 given. n 0 vn, Σn, pn and pf n being calculated. We compute ( vn+1, Σn+1) solution of

          

3 v n+1 − 4v n + v n−1 2δt + 2v n·∇v n − v n−1·∇v n−1 − ∆ v n+1 ρn+1 + ∇pn ρn+1 = div Σn+1 ρn+1 − ey, Σn+1 = Pqn(Σn+1 + r D v n+1 + ε(Σn − Σn+1)),

• v n+1
• ∂Ω

= 0.

֒ → k = 0 Σn,0 = Σn, and (vn, pn, pn

f ) given.

֒ → k 0 Σn,k known, we first compute vn,k solution of a Laplace-type problem then we project the stress tensor to obtain Σn,k+1:

          

3 v n,k − 4v n + v n−1 2δt + 2v n·∇v n − v n−1·∇v n−1 − ∆ v n,k ρn+1 + ∇pn ρn+1 = div Σn,k ρn+1 − ey,

• v n,k
• ∂Ω

= 0, Σn,k+1 = Pqn(Σn,k + r D v n,k + ε(Σn − Σn,k)).

Numerical scheme (without density)

n = 0 v0, Σ0, p0 and pf 0 given. n 0 vn, Σn, pn and pf n being calculated. We compute ( vn+1, Σn+1) solution of

          

3 v n+1 − 4v n + v n−1 2δt + 2v n·∇v n − v n−1·∇v n−1 − ∆ v n+1 ρn+1 + ∇pn ρn+1 = div Σn+1 ρn+1 − ey, Σn+1 = Pqn(Σn+1 + r D v n+1 + ε(Σn − Σn+1)),

• v n+1
• ∂Ω

= 0.

֒ → k = 0 Σn,0 = Σn, and (vn, pn, pn

f ) given.

֒ → k 0 Σn,k known, we first compute vn,k solution of a Laplace-type problem then we project the stress tensor to obtain Σn,k+1:

          

3 v n,k − 4v n + v n−1 2δt + 2v n·∇v n − v n−1·∇v n−1 − ∆ v n,k ρn+1 + ∇pn ρn+1 = div Σn,k ρn+1 − ey,

• v n,k
• ∂Ω

= 0, Σn,k+1 = Pqn(Σn,k + r D v n,k + ε(Σn − Σn,k)).

Experimental conditions

14 cm 20 cm ρ = 1550 kg.m−3 µ = 1 Pa.s δ = 34 degrees density: 1 kg.m−3 viscosity: 1 Pa.s granular bed “air”

Figure: Experimental setup

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 14 / 18

Numerical simulation

We compute the collapse with different diffusion coefficient κ. With κ = 1, it happens nothing. κ = 0.2 κ = 0.1

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 15 / 18

Diffusion coefficient = 0.2

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 15 / 18

Diffusion coefficient = 0.1

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 15 / 18

3

Perspectives

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 16 / 18

Perspectives

To compare the results given by this numerical scheme for the non-fluidised case:

◮ Lower yield stress,. . . ◮ Lower friction coefficient.. . .

To compare the results given by this numerical scheme for the fluidised case:

◮ Adapt the viscosity value. Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 17 / 18

Thank you!

Jordane Mathé (Maths-LMV) Fluidized granular flows june 2015 18 / 18