Two-phase flows with granular stress T. Gallou et, P. Helluy, J.-M. - PowerPoint PPT Presentation

Two-phase flows with granular stress T. Gallou¨ et, P. Helluy, J.-M. H´ erard, J. Nussbaum LATP Marseille IRMA Strasbourg EDF Chatou ISL Saint-Louis January 24, 2008

Context ◮ flow of weakly compressible grains (powder, sand, etc .) inside a compressible gas; ◮ averaged model (the solid phase is represented by an equivalent continuous media); ◮ 2 densities, 2 velocities, 2 pressures, 1 volume fraction; ◮ relaxation approach to return to a 1 pressure model (classical: cf . bibliography); ◮ novelty: granular stress treated in a rigorous way; ◮ stable approximation; ◮ application.

Hyperbolicity and stability Hyperbolicity Numerical viscosity A general granular flow model Two-pressure model Entropy Hyperbolicity Granular stress Reduction: one pressure model Numerical approximation: splitting method Convection step Relaxation step Numerical application Combustion chamber Conclusion Biblio

Hyperbolicity and stability

Hyperbolicity Consider w ( x , t ) ∈ R 2 solution of w t + Aw x = 0 , � 0 � ε A = , ε = ± 1 . 1 0 � ∞ −∞ e − ix ξ w ( x , t ) dx . Space Fourier transform � w ( ξ, t ) := w ( ξ, t ) = e − i ξ tA � w ( ξ, 0) . (1) � ◮ Hyperbolic case ( ε = 1): the L 2 norm of w ( ., t ) is constant ◮ Elliptic case ( ε = − 1): the frequency ξ is amplified by a factor e | ξ | t (unstable)

Numerical viscosity Classical approximations introduce a ”numerical viscosity”, which can be modelled by w t + Aw x − hsw xx = 0 . h is the size of the cells, s = ρ ( A ) the spectral radius of A . The amplification is now w ( ξ, t ) = e − i ξ tA − hs ξ 2 t � � w ( ξ, 0) . (2) But in the elliptic case, the limit system when h → 0 is still unstable ! Problem: many models in the two-phase flow community are non-hyperbolic...

A general granular flow model

◮ flow of compressible grains (powder, sand, etc.) inside a compressible gas; ◮ averaged model; ◮ 2 densities, 2 velocities, 2 pressures, 1 volume fraction; ◮ relaxation approach to return to a 1 pressure model (classical: cf. bibliography); ◮ novelty: ”rigorous” granular stress (tramway); ◮ stable approximation; ◮ application.

Two-pressure model A gaz phase k = 1, a solid (powder) phase k = 2 7 unknowns: partial densities ρ k , velocities u k , internal energies e k , gas volume fraction α 1 . Pressure law: p k = p k ( ρ k , e k ) = ( γ k − 1) ρ k e k − γ k π k , γ k > 1 E k = e k + u 2 k Other definitions: m k = α k ρ k α 2 = 1 − α 1 2 The balance of mass, momentum and energy reads m k , t + ( m k u k ) x = 0 , ( m k u k ) t + ( m k u 2 k + α k p k ) x − p 1 α k , x = 0 , ( m k E k ) t + (( m k E k + α k p k ) u k ) x + p 1 α k , t = 0 , α k , t + u 2 α k , x = ± P ,

Entropy The phase entropies satisfy the following PDEs T 1 ds 1 = de 1 − p 1 d ρ 1 ρ 2 1 T 2 ds 2 = de 2 − p 2 d ρ 2 − Θ d α 2 ρ 2 2 After some computations, we find the following entropy dissipation equation � � m k u k s k ) x = P ( m k s k ) t + ( ( p 1 + m 2 Θ − p 2 ) T 2 Natural choice to ensure positive dissipation P = 1 ε ( p 1 + m 2 Θ − p 2 ) , ε → 0 + . R := m 2 Θ is called the granular stress .

Hyperbolicity Let Y = ( α 1 , ρ 1 , u 1 , s 1 , ρ 2 , u 2 , s 2 ) T . In this set of variables the system becomes Y t + B ( Y ) Y x = S ( P ) ,   u 2 ρ 1 ( u 1 − u 2 )   ρ 1 u 1   α 1   � c 2 p 1 , s 1 u 1  1  γ k ( p k + π k )   ρ 1 ρ 1 B ( Y ) = c k =   u 1   ρ k  ρ 2  u 2   c 2 p 2 , s 2  p 1 − p 2  2 u 2 m 2 ρ 2 ρ 2 u 2 The characteristic polynomial is P ( λ ) = ( u 2 − λ ) 2 ( u 1 − λ )( u 1 − c 1 − λ )( u 1 + c 1 − λ )( u 2 − c 2 − λ )( u 2 + c 2 − λ )

Granular stress How to choose the granular stress R = m 2 Θ ? Θ = Θ( α 2 ) ⇒ Θ = 0 Thus a more general choice is necessary. Exemple: for a stiffened gas equation of state p 2 = ( γ 2 − 1) ρ 2 e 2 − γ 2 π 2 , γ 2 > 1 . We suppose Θ = Θ( ρ 2 , α 2 ). We find Θ( ρ 2 , α 2 ) = ρ γ 2 − 1 θ ( α 2 ) 2 Particular choice θ ( α 2 ) = λα γ 2 − 1 R = λ m γ 2 2 . ⇒ 2

Reduction: one pressure model When ε → 0+, formally, we end up with a standard one pressure model p 2 = p 1 + m 2 Θ (3) We can remove an equation (for example the volume fraction evolution) and we find a 6 equations system Z = ( ρ 1 , u 1 , s 1 , ρ 2 , u 2 , s 2 ) T . (4) Z t + C ( Z ) Z x = 0 . Let ∆ = α 1 α 2 + δ ( α 1 ρ 2 a 2 2 + α 2 ρ 1 c 2 1 ) , (5) and δ = α 1 − 1 /γ 2 2 . (6) λγ 2 ρ γ 2 2

Then we find The eigenvalues can be computed only numerically. We observe that when λ → 0, the system is generally not hyperbolic. We observe also that when λ → ∞ , we recover hyperbolicity.

Numerical approximation: splitting method Approximation of the one-pressure model by the more general two-pressure model. At the end of each time step, we have to return to the pressure equilibrium Relaxation approach.

Convection step Let w = ( α 1 , m 1 , m 1 u 1 , m 1 E 1 , m 2 , m 2 u 2 , m 2 E 2 ) T The system can be written w t + f ( w ) x + G ( w ) w x = Σ( P ) . In the first half step the source term is omitted. We use a standard Rusanov scheme f n i +1 / 2 − f n w ∗ i ) w n i +1 − w n i − w n i − 1 / 2 i − 1 i + G ( w n + = 0 , ∆ t ∆ x 2∆ x f n i +1 / 2 = f ( w n i , w n i +1 ) numerical conservative flux f ( a , b ) = f ( a ) + f ( b ) − s 2( b − a ) 2 For s large enough, the scheme is entropy dissipative. Typically, we take � � ρ ( f ′ ( a )) , ρ ( f ′ ( b )) s = max

Relaxation step In the second half step, we have formally to solve α k , t = ± P , m k , t = u k , t = 0 , (7) ( m k e k ) t + p 1 α k , t = 0 . Because of mass and momentum conservation we have m k = m ∗ k and u k = u ∗ k . In each cell we have to compute ( α 1 , p 1 , p 2 ) from the previous state w ∗ p 2 = p 1 + λ m γ 2 2 , m 1 e 1 + m 2 e 2 = m ∗ 1 e ∗ 1 + m ∗ 2 e ∗ 2 , ( m 1 e 1 − m ∗ 1 e ∗ 1 ) + p 1 ( α 1 − α ∗ 1 ) = 0 .

After some manipulations, we have to solve H ( α 2 ) = ( π 2 − π 1 )( α 1 + ( γ 1 − 1)( α 1 − α ∗ 1 ))( α 2 + ( γ 2 − 1)( α 2 − α ∗ 2 )) +( λα 2 m γ 2 2 − A 2 )( α 1 + ( γ 1 − 1)( α 1 − α ∗ 1 )) + A 1 ( α 2 + ( γ 2 − 1)( α 2 − α ∗ 2 )) = 0 with with A k = α ∗ k ( p ∗ k + π k ) > 0 . The solution is unique in the interval [0 , 1 − β 1 ] with β 1 = γ 1 − 1 α ∗ (8) 1 γ 1

Numerical application We have constructed an entropy dissipative approximation of a non-hyperbolic system ! What happens numerically ? Consider a simple Riemann problem in the interval [ − 1 / 2 , 1 / 2]. γ 1 = 1 . 0924 and γ 2 = 1 . 0182. We compute the solution at time t = 0 . 0008. The CFL number is 0 . 9. Data: ( L ) ( R ) ρ 1 76.45430093 57.34072568 u 1 0 0 200 × 10 5 150 × 10 5 p 1 (9) ρ 2 836.1239718 358.8982226 u 2 0 0 200 × 10 5 150 × 10 5 p 2 α 1 0 . 25 0 . 25

Figure: Void fraction, 50 cells, no granular stress .

Figure: Void fraction, 1000 cells, no granular stress .

Figure: Void fraction, 10000 cells, no granular stress .

Figure: Void fraction, 100000 cells, no granular stress . Linearly unstable but non-linearly stable...

Combustion chamber We consider now a simplified gun. The right boundary of the computational domain is moving. We activate the granular stress and other source terms (chemical reaction and drag), which are all entropy dissipative. The instabilities would occur on much finer grids...

Figure: Pressure evolution at the breech and the shot base during time. Comparison between the Gough and the relaxation model.

Figure: Porosity at the final time. Relaxation model with granular stress.

Figure: Velocities at the final time. Relaxation model with granular stress.

Figure: Pressures at the final time. Relaxation model with granular stress.

Figure: Density of the solid phase at the final time. Relaxation model with granular stress.

Conclusion ◮ Good generalization of the one pressure models; ◮ Rigorous entropy dissipation and maximum principle on the volume fraction; ◮ Stability for a finite relaxation time; ◮ The instability is (fortunately) preserved by the scheme for fast pressure equilibrium; ◮ The model can be used in practical configurations (the solid phase remains almost incompressible).