Assignment 1

One-dimensional Compressible Flow with Heat Addition and Friction


The flow inside a tube with heat addition or friction, and a stationary normal shock will be calculated in this exercise. A 1D (one dimensional) compressible flow solver will be used for the numerical simulations, 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 1-3 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_problem_3_4.m creates initial solution for problem 3.4

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 - 1D flow with friction (A)

Problem 3.12 (1 atm = 101325 Pa). To solve this problem numerically you need to calculate and change the following variables in the main script:

  • Radius of the tube (r)
  • Length of the tube (X2)
  • Friction loss coefficient (FFRIC)
  • Heat flux/unit area (QWALL)
  • Initial solution (ROINIT,UINIT,PINIT). Choose for instance the condition at the inlet of the tube.
  • Boundary conditions (BNDC.T0LB,BNDC.P0LB,BNDC.PRB) if boundary conditions of type 1 at inlet (LB) and 2 at outlet (RB) is used.

Task 2 - 1D flow with friction (B)

Problem 3.13 (1 ft = 0.3048 m, 520˚R = 288.889˚K). Try to figure out which of the previous changes that needs to be changed again. In addition to this, you need to change the boundary conditions to supersonic, since we have supersonic flow.

Task 3 - 1D flow with heat addition

Problem 3.8 (assume circular duct with length = 1m, and diameter = 1cm). Change what needs to be changed in the main script. You might need to add some extra numerical dissipation in the supersonic case (a), i.e., set DAMP.FCNTP to 0.5, and maybe use the non-reflective boundary condition at the outlet. If the computation crash even if the boundary conditions are correct, you can experiment with different initial conditions.

Task 4 - Normal shock

Problem 3.4 Here you will need to initialize the shock from the beginning and then see what happens when the simulation starts. Use the provided hint and change the variables (RO1,U1,P1,RO2,U2,P2) to your analytically calculated values. The shock should be stationary in the computational domain if everything is set correct. (assume same pipe as in task 3)

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