Bienvenue sur Code-MΦ

Un site dédié à la résolution numérique des écoulements compressibles monophasiques ou multiphasiques

Two-phase solver


PEGASE-1D computes unsteady solutions for spray equations, based on the finite volume method described in Thevand* et al..
Download the package here.

This is a FORTRAN code tested :

You can visit our page on ECOGEN website to have information about install ubuntu on Windows.
*'N. Thevand, E. Daniel, J.C. Loraud 'On high resolution schemes for solving unsteady compressible two-phase dilute viscous flows', International journal for numerical methods in fluids 31 (4), 681-702, 1999


PEGASE-1D is a free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. PEGASE-1D code is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details available at the following address: GPL license. Everyone is permitted to copy and distribute verbatim copies of this license document, but changing it is not allowed.

Using and Warning

Code-Mφ disclaims any responsibility in the obtained results with the downloaded codes.
You are encouraged to signal errors, suggest improvements to the Code-Mφ team that will try to answer.

This code is a part of the universitary program of Polytech Marseille (Master degree and diplôme d'Ingénieurs at Aix-Marseille Université). For the students, this code is necessary tool to produce a study on the interaction of a shock wave with a cloud of droplets.

Physical Model

The spray equations consits in a two-plase flow model: the carrier phase (gas phase) governed by Euler equations and the dispersed (particles or droplets) governed by pressureless Euler Equations. The two sub-systems are coupled via source terms accounting for the drag force and the heat transfer between both phases. The dispersed phase model is based on dilute flow assumptions for a monodisperse cloud of spherical particles (droplets).

Numerical scheme

This is 1D solver using regular or non-regular mesh spacing. The time integration is explicit. The global accuracy is 1rst order or higher order using MUSCL approach. The fluxes are computed thanks to an HLLC Riemann solver for the gas phase and an exact Riemann solver for the dispersed phase. This is a sequential program (run with 1 CPU).

Requirements and Installation

A Fortran compiler is required. Results can be easily plotted with Gnuplot (See

  • launch gnuplot
  • In the gnuplot terminal, enter the command line : load 'Flow.gnu'
The contains:
  • main.f90
  • This is the 1-file FORTRAN program.
  • DATAINIT INITGAZ.DAT INITPART.DAT air_gaz.inp eau_liq.inp
  • Data files contain data read by 'main.f90'.

Running PEGASE-1D

Test case: physical description

The package is ready to run and to compute the propagation of a shock wave, initially located at x=0.15 m (XODC in DATAINIT). The pressure jump is equal to 10. The unperturbated flow is at rest, the pressure is set to 1 bar and the density is equal to 1.16 kg/m3(see 'initgaz.dat' file).
The dispersed phase is made of 50 μm droplets initially at rest. The apparent density is set to 1.6 kg/m3. The temperature is set to 300.5819 (K), which is the exact gas temperature (see 'initpart.dat' file). The cloud length is 0.1 m and begins at position 0.75m (LCL and XCL in DATAINIT).
The length of the channel is equal to 1m (domain_lenght in DATAINIT).

Test case: computing data

500 timesteps are computed. Every 100 timesteps a set of results are written in file 'PLOT'. The name of this file can be changed (FILE2) in DATAINIT. This 'PLOT' file contains the flow quantities versus x in the following order (see: 'subroutine ecritplot'):
x, u, P, T, ρ, c, up, np, ρp, Tp, DIAM.
The initial solution is always written in PLOT file.

3 pressure gauges are located at x=0.35m, x=0.7d0 and x=0.99d0 (XCP1,XCP2,XCP3). In the file insta.dat (FILE1), every 1 timestep (KINSTA), are written the following quantities: x, P(XCP1), P(XCP2), P(XCP3).

The timestep obeys a stability criterion (CFL condition) and is equal to 0.5 (CFL). Nevertheless,a fixed time step can be chosen by imposing CFL=0. and a fixed value to DT0 (be sure the fixed timestep is always smaller than the CFL timestep).

The scheme order is imposed to 2 (IORDRE) with minmod limitor. Several limitors can be chosen with according the ILIMITEUR value :(0:no limitor) (1:minmod) (2:vanleeer) (3:van albada)(4:collela et woodward)(5:superbee roe)(6:superbee CHAKRAVARTHY ET OSHER)


PEGASE-1D files and structure

This a one-file code that requires no Makefile (see 'Compiling' section). The main.f90 file contains all the code.

List of modules and main features

  • MODULE precisions
  • Define the double precision type
  • MODULE dimension_probleme
  • Define the fixed array dimensions
  • MODULE courant
  • Define type for the flow quantities: cell_carrier and cell_DP
  • MODULE geometrie
  • Define type for the mesh: maillage and definition_domaine

List of function


List of subroutine (alphabetic order)

  • Compute and display the time step. The display frequency is adjusted according the total number of timesteps.
  • Euler solver for the temporal integrating of the source terms
  • Write the solution in file 'PLOT'
  • HLLC Riemann solver for the gas phase
  • 1D mesh generation
  • Flux sommation over a cell contour
  • Compute the intial condition for the gas phase. Read tags in file 'INITGAZ.DAT'
  • Compute the intial condition for the dispersed phase. Read tags in file 'INITPART.DAT'
  • Store the solution at the interafce between cells (from both Riemann solvers)
  • Boundary conditions for the gas phase (for 2nd order)
  • Boundary conditions for the dispersed phase
  • Boundary conditions for the dispersed phase (for 2nd order)
  • MUSCL8
  • Compute slope and their limitation for the gas phase
  • Compute slope and their limitation for the dispersed phase
  • Contains all the subroutines necessary for the second order
  • Projections of the vector quantities in the reference system of a cell contour
  • PRIM
  • Compute pritive quantive (i.e.: density, pressure,, velicty, temperature...) from conservative quantities (momentum, total energy...)
  • PRMN
  • Exact Riemann solver for the gas phase
  • Extrapolation of the quantities from the center to the cell contours
  • Compute physical interaction between phases (drag force, heat transfer)

Compile and run

With gfortran installed on windows 10:

  • @myfolder:> gfortran -o pegase1D main.f90
  • @myfolder:> ./pegase1D