Parareal algorithm for two phase flows simulation Katia Ait-Ameur - PowerPoint PPT Presentation

Parareal algorithm for two phase flows simulation Katia Ait-Ameur Yvon Maday (Sorbonne Université - UPMC) - Marc Tajchman (CEA) May 3, 2018 Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 1 / 23

Outline 1 Context and model Cathare numerical scheme 2 3 Numerical results Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 2 / 23

Context and model Outline 1 Context and model Cathare numerical scheme 2 3 Numerical results Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 3 / 23

Context and model Different scales of modelling Two-phase flow models (gas-liquid flows) used in the simulation of boiling in the cooling system of a nuclear power plant. Direct numerical Simulation - Meso-scale - Component scale - System scale Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 4 / 23

Context and model Motivation C ode for A nalysis of TH ermalhydraulics during A ccident and for R eactor safety E valuation Cathare essentially simulates assemblies of 1D (pipes) and 3D elements (vessels) Typical cases involve up to 10 2 or 10 3 cells with 3D elements and involve up to a million of numerical time steps Space domain decomposition method is implemented and allows a speed-up of about 4-8 using 10-12 processors Strategy of time domain decompositions, complementing the space domain decomposition, based on the parareal method Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 5 / 23

Context and model Six equation model Fine physical phenomena (description of the interfaces) are filered by the model. Flow is dominated by convection. Neglecting the viscous effects, we obtain: ∂ t ( α k ρ k ) + ∂ x ( α k ρ k u k ) = Γ k     k ) + α k ∂ x p = α k ρ k g + F int  ∂ t ( α k ρ k u k ) + ∂ x ( α k ρ k u 2  k H k + u 2 H k + u 2  � � �� � � ��  k k = α k ∂ t p + α k ρ k u k g + Q H  ∂ t α k ρ k + ∂ x α k ρ k u k   k 2 2 Main unknowns: ( p , α 1 , u k , H k ) , with α 1 + α 2 = 1 and ρ k are computed thanks to equations of state: ρ k = ρ k ( p , H k ) Γ k , Q H k : mass and energy transfers between phases F int : interfacial forces k Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 6 / 23

Context and model Closure laws Tabulated equation of state with polynomial interpolation (IAPWS) Closure laws in the momentum equations: Well posedness of the system M. Ndjinga, A. Kumbaro, F. De Vuyst, P . Laurent-Gengoux, Influence of Interfacial Forces on the Hyperbolicity of the Two-Fluid Model Interfacial friction in Cathare depends on the flow regime (bubbly, annular, dispersed,..) and on the geometry: τ i = f ( α 1 , σ, ρ 1 , ρ 2 , µ 1 , µ 2 , D h )( u 1 − u 2 ) 2 Damping term to avoid the increase of the relative velocity u r = u 1 − u 2 . Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 7 / 23

Cathare numerical scheme Outline 1 Context and model Cathare numerical scheme 2 3 Numerical results Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 8 / 23

Cathare numerical scheme Cathare scheme Simulate the components of a reactor thanks to a semi-heuristic approximation of the six-equation model Staggered mesh with scalar variables ( p , H k , α ) at cell centers normal vector ( u k ) at edges Fully implicit numerical scheme Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 9 / 23

Cathare numerical scheme Time discretisation Neglect mass and energy transfers between phases ( α k ρ k ) n + 1 − ( α k ρ k ) n  + ∂ x ( α k ρ k u k ) n + 1 = 0   ∆ t    ( α k ρ k u k ) n + 1 − ( α k ρ k u k ) n   k ) n + 1 + α n + 1 ∂ x p n + 1 = ( α k ρ k ) n + 1 g + F n , n + 1  + ∂ x ( α k ρ k u 2    k k ∆ t   � � n , n + 1 � n − 1 , n � H k + u 2 H k + u 2 1 � � ( α k ρ k ) n + 1 k − ( α k ρ k ) n k   ∆ t 2 2      �� n + 1 p n + 1 − p n  H k + u 2 � �  = α n + 1 + ( α k ρ k u k ) n + 1 g  k + ∂ x α k ρ k u k   k 2 ∆ t  Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 10 / 23

Cathare numerical scheme Space discretisation Upwind scheme to express cell centered unknowns at edges in mass and energy equations. ( α k ) i − 10 − 5 , ( u k ) i + 1 / 2 > 0 � ( α U k ) i + 1 / 2 = ( α k ) i + 1 − 10 − 5 , ( u k ) i + 1 / 2 < 0 � ( ρ k ) i , ( α U k ) i + 1 / 2 ( u k ) i + 1 / 2 > 0 ( ρ U k ) i + 1 / 2 = ( ρ k ) i + 1 , ( α U k ) i + 1 / 2 ( u k ) i + 1 / 2 < 0 The convection term: � ( u k ) i + 1 / 2 (( u k ) i + 1 / 2 − ( u k ) i − 1 / 2 ) , ( u k ) i + 1 / 2 > 0 ( u k ∂ x u k ) i + 1 / 2 = ( u k ) i + 1 / 2 (( u k ) i + 3 / 2 − ( u k ) i + 1 / 2 ) , ( u k ) i + 1 / 2 < 0 Semi-heuristic approach: if ( u k ) i + 1 / 2 ≤ ( u k ) i − 1 / 2 and ( α k ) i ≤ 10 − 3 then: � ( u k ) i + 1 / 2 − C 1 ( α )( u k ) i − 1 / 2 + C 2 ( α )( u k ) i + 1 / 2 � ( u k ∂ x u k ) i + 1 / 2 = ( u k ) i + 1 / 2 C 1 ( α ) + C 2 ( α ) Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 11 / 23

Cathare numerical scheme Non linear solver Newton scheme: The semi-discretised problem: U n + 1 − U n + A ( U n + 1 , U n ) = S ( U n ) ∆ t ∆ U k + 1 + J ( U k + 1 , U k )∆ U k + 1 = S ( U n , U k ) , where: ∆ U k + 1 = U k + 1 − U k ∆ t and : U k + 1 = ( P k + 1 , α k + 1 , H k + 1 , H k + 1 , u k + 1 , u k + 1 ) V l v l v In Cathare: By Gauss elimination, obtain a system with pressure increment only. Solve the problem in pressure with a direct linear solver (LAPACK BLAS) In MiniCathare (Cathare restricted to 1 test case): solve the complete linear system with an iterative linear solver (PETSC library) Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 12 / 23

Cathare numerical scheme Characteristics of Cathare scheme Accuracy of the scheme for nearly incompressible flows For single phase flows: Riemann solvers have poor precision in the incompressible limit S. Dellacherie, Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number Staggered schemes enjoy good precision at the incompressible limit Two phase flows: Special treatment of the vanishing phase M. Ndjinga, T. P . K. Nguyen and C. Chalons, A 2x2 hyperbolic system modelling incompressible two phase flows: theory and numerics Countercurrent flows Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 13 / 23

Numerical results Outline 1 Context and model Cathare numerical scheme 2 3 Numerical results Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 14 / 23

Numerical results Oscillating manometer G.F . Hewitt, J.M. Delhaye, N. Zuber, Multiphase science and technology Ability of a scheme to preserve system mass and to retain the gas-liquid interface Flow regime: separated phases Interfacial friction term to handle the vanish- ing phase and adaptation of the convection term Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 15 / 23

Numerical results Stopping criteria Initial condition: P = 10 5 , h l = 4 . 17 × 10 5 , h v = 2 . 68 × 10 6 , u v = u l = − 2 . 1 and � 1 − 10 − 5 , in the upper half α v = 10 − 5 , elsewhere Time interval : [0,20] Order of convergence in time of the Cathare scheme: Reference solution: 220 cells and δ t = 10 − 5 Error norm: max n || U n − U n ref || L 2 where U n = P n , α n v , h n v , h n l , u n v , u n � � l max n || U n ref || L 2 Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 16 / 23

Numerical results Order of convergence in time Error on 110 cells Relative error : L 2 in space, L ∞ in time Error on 220 cells √ C ∆ T 0 . 4 0 . 3 0 . 2 0 . 1 0 . 0002 0 . 0004 0 . 0006 0 . 0008 Time step D. Bouche , J.-M. Ghidaglia , F. Pascal, Error estimate and the geometric corrector for the upwind finite volume method applied to the linear advection equation Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 17 / 23

Numerical results Parareal for hyperbolic equations J.-L. Lions , Y. Maday , G. Turinici , Résolution par un schéma en temps "pararéel" Initialisation: k = 0 , U 0 n + 1 = G ( T n , T n + 1 , U 0 n ) sequential Parareal iteration k: ( U k n ) N n = 0 known. (k.1) Compute fine solution on each ] T n , T n + 1 [ : F ( T n , T n + 1 , U k n ) in parallel (k.2) Prediction coarse step: G ( T n , T n + 1 , U k + 1 ) sequential n (k.3) Correction step: U k + 1 n + 1 = G ( T n , T n + 1 , U k + 1 ) + F ( T n , T n + 1 , U k n ) − G ( T n , T n + 1 , U k n ) n Convergence properties [Gander, Vandewalle, 2007] and stability analysis [Maday, Ronquist, Staff, 2005] Dependence on the regularity of the initial condition and solution [Bal, 2005] and [Dai, Maday, 2013] Correction procedure re-using previously computed information based on a projection on a Krylov subspace [Gander, Petcu, 2008] or on a reduced basis [Chen, Hesthaven, Zhu, 2014] . Coarsening in space [Ruprecht, 2014] and [Lunet, 2017] Katia Ait-Ameur (CEA - UPMC) 7th PinT Workshop May 3, 2018 18 / 23