Assignment 3

Quasi-One-Dimensional Flow


Flow inside a nozzle, and unsteady wave motion will be studied in this exercise. A quasi 1D compressible flow solver will be used for the numerical simulations (same solver as in Assignment 1) and the results should be compared with theory. A Matlab code needed for this assignment g1deul_MAIN.m and associated Matlab functions are available in the code archive below.

OBS Read chapters 5-7 in the course book before starting with the exercise!

When you have done the assignment, you should write a short report and hand in for assessment of this course element (instructions at the bottom of this page).
Code Archive
g1deul_MAIN.m main script
mtr.m function that computes inverse of cell volume
stepexpl.m function performing explicit time step
auxvr.m function that calculates auxiliary variables
hint_harmonic_wave.m creates initial solution for harmonic wave problem
hint_problem_5_11_nozzle.m Q1D discretization for convergent-divergent nozzle

G1DEUL manual

The code G1DEUL is based on FVM (finite volume method) and can be used to study one-dimensional stationary or transient compressible flow phenomenon. Basic information of how to run the code G1DEUL in MatLab is given below.

Only the main script g1deul_MAIN.m needs to be edited in this exercise. The other files should not be changed. When the main script is executed all the initializations of variables needed for the solver is done and the main time loop starts. The solution is plotted in a window during the calculation, but this is not done every time step since it takes a lot of extra time. The plotting interval can be changed by pressing the buttons “+10” or “-10” in the plot window. The window contains plots of flow variables \(\rho\) (density), \(u\) (velocity), \(c\) (speed of sound), \(p\) (static pressure), \(T\) (static temperature), \(M\) (Mach number), \(p_o\) (total pressure), \(T_o\) (total temperature). It also contains plots of the residual \(R_{max}\). The cross section area of the tube (\(SX\)) is plotted in a separate window.

Matlab call-back window

The CFL-number is for compressible flow defined as:

$$CFL=\frac{\Delta t(c+|u|)}{\Delta x}$$

where \(\Delta t\) is the time step size, \(c\) is the speed of sound, \(u\) is the fluid velocity, and \(\Delta x\) is the local cell size. Since the solver uses explicit time marching, it is important that the CFL-number is below unity, otherwise the computation may diverge. The time step size is calculated in the main loop from a specified CFL-number.

The variables in the main script are organized into structs. After running the main script, the value of a variable can be printed in the "Command Window" by typing the name of the struct followed by a dot and then the name of the variable;


ans =


To see all the variables in a structure, type only the name of the structure;

>> MSH

MSH = 

      NI: 101
      XM: [101x1 double]
      SX: [101x1 double]
    SCFL: [101x1 double]
    VOLI: [100x1 double]
     XCP: [100x1 double]

The structure MSH contains information about the mesh; NI = number of mesh points, XM = vector of mesh points (x-coordinates), SX = vector of cross section areas, SCFL = vector of circumferential lengths, VOLI = vector of inverse of cell volumes, and finally XCP = vector of cell center x-coordinates;

mesh nodes and cell centers

The variable NI and the vectors XM and SX should be specified. Then the vectors XCP, VOLI, and SCFL is calculated (assuming a circular pipe for the last one). The vectors SX, and SCFL corresponds to the mesh points vector XM, and the vector VOLI corresponds to the cell centers vector XCP.

The other structures that needs to be initialized before the solver can be used are: GAS, DIFF, SOLV, BNDC, COEF, DAMP, and TSTP

GAS = 

    GCNT: 1.4000
    RGAS: 287

The structure GAS contains the properties of the gas, where GCNT = \(\gamma\) (gamma), and RGAS = \(R\).


    FFRIC: 0
    QWALL: 213934

The structure DIFF contains the wall friction and heat addition at the walls of the tube; FFRIC = wall friction (dimensionless), and QWALL = heat flux per unit area (W/m2). Note that the definition of the heat flux here is not the same as in the book, where it is defined as heat added per unit mass (J/kg). The following formula can be used to convert between the different definitions


where \(\dot{q}_{wall}\) is the heat flux per unit area (W/m2), \(q_{mass}\) is the heat added per unit mass (J/kg), \(\dot{m}\) is the mass flow (kg/s), \(l\) is the total length of the tube (m), and \(b\) is the circumferential length of the tube (m).


      RO: [100x1 double]
      RU: [100x1 double]
      ET: [100x1 double]
       U: [100x1 double]
       P: [100x1 double]
       C: [100x1 double]
    TIME: 0.0819

The structure SOLV contains the flow variables used by the solver, at each cell center. An initial solution of RO (\(\rho\); density), RU (\(\rho u\); density times velocity), and ET (\(\rho e_o\); density times total energy) has to be created and the TIME variable has to be initialized before the calculation can start. The other variables U (\(u\); velocity), P (\(p\); pressure), and C (\(c\); speed of sound) are calculated by the function "auxvr".


    ITLB: 1
    ROLB: NaN
     ULB: NaN
     PLB: NaN
    T0LB: 275.2000
    P0LB: 1.0414e+005
    ITRB: 2
    RORB: NaN
     URB: NaN
     PRB: 7.2733e+004
    T0RB: NaN
    P0RB: NaN

The BNDC structure contains information about the left and right boundaries. There are several choices of boundary conditions, and the variable ITLB defines what boundary condition there is on the left side, and ITRB defines the right side boundary condition. In this example the left boundary is of type one, i.e. subsonic inflow, and then the total temperature (T0LB) and the total pressure (P0LB) are set to 275.2 K and 104136.4 Pa respectively. The right side boundary is set to type two, i.e. subsonic outflow, and then the static pressure (PRB) is set to 72733.4 Pa. The other variables can be set to anything (here NaN; not a number) and it will not affect the solution. The available boundary conditions are;

IT	RO	U	P	T0	P0	boundary condition
0	used	used	used	-	-	non-reflective
1	-	-	-	used	used	subsonic inflow
2	-	-	used	-	-	subsonic outflow
3	used	used	used	-	-	supersonic inflow
4	-	-	-	-	-	supersonic outflow
5	-	-	-	-	-	solid wall
10	-	-	-	-	-	periodic

Boundaries of type 0-4 will be used in this exercise. As the names imply, different boundary conditions should be used if the flow is subsonic or supersonic at the boundary. The non-reflective option can be used with both subsonic and supersonic flow at the boundary.


    COEF1: -0.1667
    COEF2: 0.8333
    COEF3: 0.3333
    COEF4: 0

The structure COEF contains the coefficients for the numerical flux scheme. The default scheme is a standard third order upwind scheme.


    FCNTP: 0
    FCNTC: 0

The structure DAMP contains coefficients for extra pressure- and density-controlled numerical dissipation. If compression shocks are present in the solution the FCNTP variable can be set to 0.5 to smooth out oscillations around the discontinuity. The variable FCNTC can be set to 2 if strong contact discontinuities are present.


      CFL: 0.9000
    NSTEP: 10000
       DT: 8.2293e-06

The TSTP structure contains the CFL-number, the number of time steps (NSTEP), and the current time step size (DT).

All of the following tasks should be solved both numerically and analytically (if possible). Note that in order to be able to apply the correct boundary conditions in the numerical simulations, you might need to solve the problems analytically first! Report graphs of \(\rho\), \(u\), \(p\), \(T\), \(M\), \(p_o\), \(T_o\), from the numerical simulations and explain what you see. You should also explain what boundary condition that is used for each numerical simulation (type and values), along with other parameters such as time step size, number of time steps etc.

Task 1 - Quasi-1D simulation of the flow in a convergent-divergent nozzle

Problem 5.11

Run different variants of Problem 5.11 in the course book (assume \(T_o\) = 288°K for the numerical simulation). Use the provided hint (code archive above) when generating the nozzle geometry. Use appropriate boundary conditions for the three different cases. Compare your simulation results with your calculations.

  1. Setup Problem 5.11 in the quasi-1D compressible flow solver and run it. How do the results compare with your calculations? Compare e.g. A/A*.
  2. Calculate the pressure which results in Mach=1 at the throat but subsonic in the entire divergent part of the nozzle. Run the case in the quasi-1D flow solver.
  3. Calculate the pressure below which, results in supersonic flow throughout the entire divergent part of the nozzle. Run the case in the quasi-1D flow solver.

Task 2 - shock tube

Problem 7.10

Use a constant area tube with 200 cells and solid walls on both boundaries and make a piecewise constant initial solution (as you did in assignment 1). Assume \(p_4 = 500 000\) Pa for the numerical simulation.

Solve the problem analytically and compare with the results obtained in Matlab.

You should calculate and compare

  • The shock speed
  • The head and tail speed of the expansion wave
  • The velocity of the gas behind the shock
  • The velocity of the reflected shock

One way of finding the shock, head, and tail velocities is to pause the computation after a few time steps and note the position and simulation time. The speed can be estimated from the distance traveled from the center of the tube. Use (DAMP.FCNTP=0.5, DAMP.FCNTC=2) in this case since there will be both strong shocks and contact discontinuities.

If you think that the animation of the shock movement in the figure is too fast, you can change the variable button.plotstep to e.g. 5 (default setting is 20). Since plotting is highly time consuming compared to the calculation time for this program, this change will slow down the animation.

Task 3 - traveling acoustic wave

You should make a simulation of a harmonic acoustic wave in a constant area tube. Use periodic b.c. at both sides and initialize a harmonic wave (using the hint on the course web page). Let the wavelength, \(\lambda\), be the length of the tube. The tone should be a standard A440 (440 Hz), and let it travel for 10 m. At approximately what sound pressure level does nonlinear effects become significant? Can this be estimated analytically, and if so how?


The relation between the wave length \(\lambda\) (m), the the speed of sound in the undisturbed fluid \(a_{\infty}\) (m/s), and the frequency \(f\) (Hz) is as follows:

The sound pressure level \(L_p\) (dB) is calculated according to:


where \(p_{rms}\) is the root-mean-square of the pressure signal, and \(p_{ref}\) is a reference pressure (here we used the value \(p_{ref}=2 \times 10^{-5}\) Pa).

For harmonic waves the equation above can be rewritten as:

$$L_p=20\log_{10}\left(\frac{\Delta p}{p_{ref}\sqrt{2}}\right)$$

where \(\Delta p\) is the amplitude of the pressure harmonic wave (Pa).


  • Stop the calculation after the wave has traveled 10m by adding the following code lines at the end of the main for-loop.
	 if SOLV.TIME > 10/Cinf
  • To compare the profile of the sound wave for different sound pressure levels, save the velocity field (SOLV.U) and the coordinate vector (MSH.XCP) for a number of different values. Normalize the velocity with the maximum value for each case and gather the profiles in one figure.

Assignment Report

In order to complete the assignment, a report should be written and handed in for assessment. Add a section in the report for each of the tasks above including:

  1. Given data
  2. Analytical solution
    • Hand-written calculations can be included as images
  3. CFLOW solver settings
    1. Geometry settings
    2. Mesh properties
    3. Numerics
      • numerical scheme
      • artificial damping (if applicable)
      • time integration
    4. Run parameters
    5. Gas properties
    6. Initial flow field
    7. Boundary conditions
      • selected boundary condition type and specified values for both left and right boundary
    8. Sources
  4. Solver convergence
    • Residual plot
    • Do you consider the solution to be converged enough? Justify your answer
    • Where there any problems getting the solver to converge for this specific case. If so, speculate on the cause for this behavior.
  5. Results
    • Flow field plots (Mach number, Temperature, Pressure)
    • Additional flow field plots if it shows something of interest (for example total pressure and/or total temperature)
    • Comment on what you see (specific flow features that can be observed)
    • Do the results make sense?
  6. Discussion
    • Are there any significant deviations from the analytical/theoretical solution? If so, what is the reason for the deviations? (Might be difficult to answer exactly but a speculation will do)
    • Are there any indications of numerical errors in the solution?
      • wiggles in flow field data
      • changes in total temperature for cases that should be adiabatic
      • changes in total pressure for cases that should be isentropic
      • other indications of numerical errors
      If you see anything in your solutions that you think is an indication of numerical errors, try to explain why those errors appears.

When done, submit the report in Canvas