Visualizing 2D Shock Wave Propagation Inside A Shock Tube Channel Using Converge CFD Software & Post-processing It In ParaView

Objective:- The experiment is focused on setting up a shock tube simulation using a transient flow solver to visualize shock waves generated due to the pressure difference created between the two regions when the diaphragm inside the tube is ruptured. Such an experiment is primarily useful in the study of chemical kinematics where a part of the tube is filled with fuel with the other section filled with air & due to these shockwaves there would be auto ignition & from these data, we can determine the time parameters involved for the fuel & air to combust. This is known as 'ignition delay'. For this particular project, we will not be creating a diaphragm as this would increase the complexity of the challenge instead we would be making use of "events" in Converge to specify whether the regions (high pressure & low pressure) are connected or not in a particular time step. Our selection of species for the project is nitrogen (N2) & oxygen (O2). We will be setting up a simple 2D simulation, which will only involve mixing between two species and no chemical reaction. Through this project we will also analyze the solution as a function of time to see how it behaves as time progresses. We will also run multiple cases for different grid resolutions and analyzed how the number of cells considered affects the solution. Before setting the case setup we will be also making use of boundary flagging techniques for assigning each side of the geometry to a particular boundary & these boundaries will later be assigned to a volumetric region. Next, we set up the simulation parameters using a transient solver. Once the setup is complete, we export these inputs files generated by Converge CFD & run it using CYGWIN to generate 3D post output files. These 3D output files are first converted into either Paraview vtk inline binary format or EnSite format for Paraview to read & generate the flow field simulation. 

                                                                                                                                                                                                                                                                     

                                                                                                 Fig 1.1:- Sketch of a shock tube at initial condition.

Physical Problem Description                                                                                                                                                                                                                                                                  

A shock tube consists of a channel tube in which a gas at high pressure, called the driver gas, is initially separated from gas at a lower pressure, called the test gas or the driven gas, by a diaphragm. The high and low-pressure regions can also be referred to as compression and expansion pressure chambers respectively. The initial state is defined by the values for density, pressure, and velocity, as shown in Figure 2.1. All the viscous effects are negligible along the tube sides; it is also assumed that there is no motion at the beginning of the experiment. The initial state at time t = 0 consists of two constant states 1 and 2 with p2 > p1, t2 = t1, and v1 = v2 = 0. Just after the membrane is removed (for our case we will be using the 'event' feature to control the opening of the chamber membrane) an incident shock wave generates and propagates through the test gas raising its temperature and pressure. As the shock wave moves through the test gas, a rarefaction wave moves back into the high-pressure gas at the speed of sound and reflects off the driver's end wall once it reaches it. The test gas and the driver gas make contact at the “contact surface”, which moves along the channel tube behind the shock front. The incident shock wave arrives at the end wall and reflects back raising the temperature and pressure of the shocked gas in the test region. As such, the energy release and chemical reaction start resulting in an ignition event. The end of the experiment is dictated by the arrival of the expansion fan to the driven tube end wall or by the interaction of the reflected shock with the contact surface. Typical test times in the shock-tube experiments are on the order of few milliseconds. Fig 2.2 shows the ideal movement of the shock front, the contact surface, the rarefaction wave, and the reflected shock wave in a distance-time diagram.  

                                                           

    

                                                                                          Fig 2.1:- Sketch of a shock tube at the initial state with discontinuities.

                                                       

 

                                                                 Fig 2.2:- Plot shows different waves that are formed in the shock tube once the diaphragm is ruptured.

Shock Tube Geometry Setup                                                                                                                                                                                                                                                                  

The pre-modeled channel tube was designed using Solidwork’s & imported into Converge Studio as an STL (stereolithography or standard triangle language) file. An important fundamental concept of Converge software is that any geometry designed in converge or any geometry exported to converge will be made up entirely of triangles. After importing the STL file of the duct into Converge we convert the dimensions. Since, the geometry was created in solid works, the dimensions of the geometry are expressed in millimeters or ‘mm’ but the Converge software package uses S.I unit systems i.e. ‘meters’ & thereby assumes that the geometry imported is in meters. To fix this we use the geometry editing tools available in converge, we ‘Transform’ the entire geometry by a uniform scale factor of 0.001 as shown in Fig 3.1.

Once, the geometry is scaled down we run a ‘Diagnostic’ test to check for any anomaly like intersection errors, nonmanifold problems, open edges, etc within the geometry contour. For our case we get a number of non-manifold edges. A shock tube consists of 2 regions, a high-pressure region, and a low-pressure region and a diaphragm separating the two as stated earlier. The diaphragm is also present in our geometry, which causes the non-manifold edges error, so it is deleted as in our case we will be using a virtual diaphragm using the 'events' option provided by Converge. If we run a second diagnostic test after deleting the diaphragm we observe there are no errors, which is denoted by “green checks”, then we proceed to check for the orientations of the “Normal’s” in the geometry.

Every geometry will have a normal vector which is perpendicular to the geometry, for this problem we use the “Normal Toggle” option to check for the direction of the normal’s. If the normal’s are pointing outside of the geometry, then it’s essential to transform these normals to point inside the geometry, where the fluid flow will occur as shown in Fig 3.1. 

In Converge, all the surface information is stored in the form of 3 entities vertices, edges & triangles. Therefore, any geometry created or exported from other CAD software (like in our case) into converge is converted into triangles & that is the fundamental entity in Converge, where every part of the geometry is assumed to be a triangle. The process of grouping these triangles into boundaries is known as boundary flagging. Each boundary is assigned with a distinct ID as shown in Fig 3.2. For our geometry we have 2 distinct surfaces High_pressure & Low_pressure region.

 

     

                                                                                    Fig 3.1:- Geometry dimensions, normal orientation & diagnostic test with no errors.

                                                  

                                                                                                                                 Fig 3.1:- Boundary Flagging.

Setting Up The Flow Physics For The Computation Model                                                                                                                                                                                                                  

The simulations were carried out with the Converge CFD Software. In Converge after editing the geometry, the next stage of the process is to set up the simulation parameters. This step involves specifying certain properties to capture the physics of the problem perfectly, such as materials used, simulation time parameters, solver parameters, initial & boundary conditions to solve the Navier Stokes equation which are PDE’s, defining body forces, type of fluid flow, species involved & grid size of the mesh.

We start setting up the time-based general flow simulation by first defining the Materials involved in the simulation. Since, we are simulation 2 different fluid flow inside a channel we choose ‘gas simulation’ and our choice of species are nitrogen (N2) & oxygen (O2), & we choose our pre-defined mixture as air. For gas simulation parameters we will be using the Redlich-Kwong equation of state with a critical temperature of 133K & critical pressure of 3770000pa. For global transport parameters we use the default values of ‘turbulent Prandtl number’ and ‘Schmidt number’ which is 0.90 & 0.78 respectively.

The next parameter for our case setup is the Simulation Parameters, which include the ‘run parameters’, ‘simulation time parameters’ & the ‘solver parameters’. For the shock tube problem, we will be running compressible gas flow using transient solver at full hydrodynamics simulation mode, since our geometry is simple & has no moving part & we also require to solve the NS equations. In addition to this we will be using ‘density-based’ PISO Navier strokes solver. As for simulation time parameters we will be running the simulation for an end time of 0.003 with an initial & minimum time step as (1e-09).

As earlier, after designing the geometry, we grouped the triangles and assigned them to a particular boundary ID. Similarly, in Boundary Condition, we group the 2 boundaries & assign each of them to a volumetric region. For our case we add two volumetric regions (region 0 & region 1) from Initial Conditions & Events & assign the volumetric region Region 0 to High_pressure boundary while we assign the volumetric region Region 1 to Low_pressure boundary. These volumetric regions are required to set up the initial conditions for solving PDEs. Initial conditions are assigned at the volume whereas boundary conditions are given at boundaries. The High_pressure regions are stationary slip walls with an initial pressure of 6000000 pa with temperature conditions specified as 300K with Nitrogen(N2) as the species present inside the region. The Low_pressure regions are stationary walls with an initial pressure of 101325 pa & temperature conditions of 300K with Oxygen (O2) as the species present inside the region.

As we have already mentioned earlier our project will not contain a physical diaphragm instead we will be creating a virtual diaphragm by using a feature called Events. The time period when the 2 regions should remain closed(at 0 Seconds) and when it should open(at 0.001 Seconds) are specified.

                                                                                         

The final parameter that we need to setup is defining the Mesh that is the Base grid size. For our case mesh grid size of dx=dy=dz=0.002 & dx=dy=dz=0.003 will be compared. In addition, we will also be using a species-based Adaptive Mesh Refinement (AMR) technique to monitor the temperature at certain locations where a shock is expected to occur. Adaptive Mesh Refinement was provided with the maximum embedding of 3 layers & a sub-grid scale of 0.001. Sub-grid scale (SGS) is a parameter that will be used by Converge for refining the mesh at a specific location. Converge monitors the curvature of the fluid property provided in the AMR and refines the mesh when the property curvature variation between the consecutive grid is more than the defined SGS value. We will be running the simulation for a specific case of SGS keeping the base mesh size constant.

Post Processing                                                                                                                                                                                                                                                                                        

The function of the Converge studio is to set up the simulation & then create several input files which are then exported to a particular folder. To run these input files we use CYGWIN, a command-line interface that reads these inputs files & solves the complex PDEs of governing equation to generate several output files. The simulation is run in parallel with 4 processors. These output files are then post-converted into 3D output files which are readable files for ParaView.

Observations & Inference for Case -2 (Base Grid of 0.003)                                                                                                                                                                                                                 

 

  A. Computational Mesh With AMR

                                                                                                                                           Initial Timestep

                                                                       

                                                               Time step = 25                                                                                                                                           Time step = 50

     

                                                                Time step = 75                                                                                                                                         Time step = 100

   

                                                             Time step = 125                                                                                                                                          Time step = 150

   

                                                               Time step = 175                                                                                                                                     Time step = 200

   

                                                                                                                                         Final Time step

                                                            

 

                                                                         

                                    

    

The total number of cells generated for a mesh grid size of dx=dy=dz=0.003 during the initial time steps is lower until 0.0011 seconds as the time steps increased & the virtual diaphragm was ruptured causing the fluid (N2) to leave the high-pressure region to the low-pressure region of the shock tube. This movement of the fluid from the high-pressure region to the low-pressure region was captured by using the species based AMR with 3 layers of embedding. In the case setup we made use of species AMR for N2, which resulted in a 3 layer refinement near the regions of the fluid N2 motion. The AMR & Events feature was enabled after 0.001 seconds, hence we start observing the high number of cell generation after 0.001 seconds after which leads to a non-linear change in the number of cell generation. The animation shows the mesh development for the entire geometry from the initial time step to the final time step. 

 

  B. Pressure Contour & Mass Fraction Variation

                                                                                                                                                Initial Time Step

                                                              

                                                   Time Step = 25                                                                                                                                                   Time Step = 50

             

                                                   Time Step = 75                                                                                                                                                  Time Step = 100

             

                                                   Time Step = 125                                                                                                                                                Time Step = 150

             

                                                   Time Step = 175                                                                                                                                               Time Step = 200

             

                                                                                                                                            Final Time Step

                                             

                                                                            

 

                                       

  

                                                                                                                                            Initial Time Step

                                     

                                                                        Time Step = 25                                                                                                                    Time Step = 50

         

                                                                       Time Step = 75                                                                                                                    Time Step = 100

        

                                                                      Time Step = 125                                                                                                                   Time Step = 150

        

                                                                       Time Step = 175                                                                                                                   Time Step = 200

       

                                                                                                                                         Final Timestep

                                        

 

                                                                   

 

At the initial time step for both Pressure Contour & Mass Fraction (N2), the red section on the left of the tube is the High_pressure region with nitrogen (N2) inside it & the blue section on the right of the tube is filled with oxygen (O2) which is the Low_pressure region. In case setup we initially specified 600000pa of pressure in the high-pressure region with 101325pa in the low-pressure region.  For the initial time steps, we can observe the high-pressure region & low-pressure region have been contained in separate sections of the tube successfully by a virtual diaphragm which indicates we have used our "Events" feature correctly. Also, we can refer to the pressure plot to verify the same, from the pressure plot we observe a constant plot until 0.0011s after which it varies in a cyclic manner as the flow furthers. The varying of the plot starts as soon as the virtual diaphragm is opened i.e. after 0.001 seconds as we have mentioned in the case set up, during this time period the diaphragm will open & we observe the mixing of the two entities inside the shock tube. We observe at time step=25 the diaphragm breaks, a shock, slip surface, and expansion waves are formed and propagate through the tube. A formation of a normal shock wave, which propagates along the length of the tube until the end of the driven section which is observed at time step=50. Then it gets reflected by the wall and hence returns back, which can be noted as the mass fraction of N2 gas retracts back, and the wave furthers towards the left end. The same is depicted in the pressure contour when the right part of the tube is observed precisely. We observe when the shock wave is created, the expansion or rarefaction waves are created at the same instant. These expansion waves move in the direction opposite to that of the shock wave with the speed of sound causing continuous or smooth changes in the state of the gas. Then they get reflected by the end wall of the driver side and cause the expansion along the length of the shock tube as they move. [Note: The expansion and other changes in the state of gas happen in the direction opposite to the direction of propagation of expansion wave]. Then again, as the shock wave, after reflection from the left end, travels through the length of the tube causing the gases moving leftward to move rightward to an extent. Meanwhile, the contact surface has moved through a distance from its initial position due to the motion of the gases and the influence of pressure. Based on the above inferences the shock tube could be divided into 4 regions as shown below: 

                                                                                 

                                                                                                                  Fig:- Shock tube shortly after diaphragm has burst.

• • Region 1: It is the region between the shock and the right end of the wall, where the gas is at rest and the shock wave has not affected it yet.
  • Region 2: It is the region between the shock wave and the contact surface, which is affected by the initial shock wave.
  • Region 3: It is the region between the contact surface and the tail of the expansion fan, which is under expansion & the gas is moving in the direction of the shock.
  • Region 4: It is the region between the head of the expansion wave & the left end of the wall, where the gas is at rest & yet to be affected by the expansion wave.  
   

C. Velocity Contour Variation Inside The Shock Tube

                                                                                                                                                  Initial Timestep

                                                      

                                                                       Time Step = 25                                                                                                                       Time Step = 50

 

                                                                       Time Step = 75                                                                                                                    Time Step = 100

 

                                                                     Time Step = 125                                                                                                                   Time Step = 150

  

                                                                    Time Step = 175                                                                                                                   Time Step = 200

  

                                                                                                                     Final Timestep

                                   

                                                                

At the initial time step for velocity Magnitude, the High_pressure region with nitrogen (N2) & the Low_pressure region with oxygen (O2) is static indicating that there is no visible movement throughout the tube crosssection as the virtual diaphragm is still closed. Once the diaphragm is opened after 0.001 seconds we start observing velocity variation throughout the shock tube cross-section. As we can observe from the animation the variation is very drastic throughout the domain. At time step=25, after the diaphragm is opened we observe the initial development of the velocity profile propagating on both sides of the tube. As the time steps increased the drastic variation in velocity is attributed to the formation of shock waves. However we observe there is no considerable change in velocity of the gases near the ends, which denotes that the gases are at rest at these locations. Throughout the entire time steps, we observe the velocity profile to be propagating only in the center of the tube.

 

D. Temperature Contour Variation Inside The Shock Tube

                                                                                                                                                  Initial Timestep

                                   

                                                                         Time Step = 25                                                                                                                   Time Step = 50

 

                                                                          Time Step = 75                                                                                                                   Time Step = 100

 

                                                                           Time Step = 125                                                                                                                   Time Step = 150

 

                                                                         Time Step = 175                                                                                                                   Time Step = 200

 

                                                                                                                                    Final Timestep

                                  

                                                                       

                                    

 

At the initial time step for temperature, in the case setup we initially specified 300 K of temperature for both the high-pressure & low-pressure regions. We can observe this from the animation & the temperature plot as in the initial time steps there is no variation & we observe a straight line. As the diaphragm is broken at 0.001 seconds we start observing a cyclic motion of temperature variation across the cross-section of the tube. There is a smooth variation in a particular section and a drastic variation in some sections. The contact surface and its motion can be identified clearly from the above contour by means of the temperature discontinuity prevailing in the contour in the left section of the tube. As the shock wave causes a drastic increase in the temperature as it passes by, the contact surface is found to have completely different projections of temperature when compared to the other regions. 

 

Observations & Inference for Case -2 (Base Grid of 0.002)                                                                                                                                                                                                                      

 

  A. Computational Mesh With AMR

                                                                                                                                           Initial Timestep

                                          

                                                               Time step = 25                                                                                                                                           Time step = 50

                      

                                                                Time step = 75                                                                                                                                         Time step = 100

                  

                                                               Time step = 125                                                                                                                                          Time step = 150

                  

                                                               Time step = 175                                                                                                                                          Time step = 200

                  

                                                                                                                                         Final Timestep

                                               

 

                                                                         

                                                                     

                                                                     

    

The total number of cells generated for case 2 > case 1 even during the initial time steps is lower until 0.0011 seconds as the time steps increased & the virtual diaphragm was ruptured causing the fluid (N2) to leave the high-pressure region to the low-pressure region of the shock tube. This movement of the fluid from the high-pressure region to the low-pressure region was captured by using the species based AMR with 3 layers of embedding. In the case setup we made use of species AMR for N2, which resulted in a 3 layer refinement near the regions of the fluid N2 motion. The AMR & Events feature was enabled after 0.001 seconds, hence we start observing the high number of cell generation after 0.001 seconds after which leads to a non-linear change in the number of cell generation. The animation shows the mesh development for the entire geometry from the initial time step to the final time step. 

 

  B. Pressure Contour & Mass Fraction Variation

                                                                                                                                                Initial Time Step

                                              

                                                   Time Step = 25                                                                                                                                                   Time Step = 50

             

                                                   Time Step = 75                                                                                                                                                  Time Step = 100

             

                                                   Time Step = 125                                                                                                                                                Time Step = 150

             

                                                   Time Step = 175                                                                                                                                               Time Step = 200

             

                                                                                                                                            Final Time Step

                                             

                                                                            

 

                                                                                                                                                                                                                         

  

                                                                                                                                            Initial Time Step

                                     

                                                                        Time Step = 25                                                                                                                    Time Step = 50

         

                                                                       Time Step = 75                                                                                                                    Time Step = 100

        

                                                                      Time Step = 125                                                                                                                   Time Step = 150

        

                                                                       Time Step = 175                                                                                                                   Time Step = 200

       

                                                                                                                                         Final Timestep

                                        

 

                                                                   

 

At the initial time step for both Pressure Contour & Mass Fraction (N2), the red section on the left of the tube is the High_pressure region with nitrogen (N2) inside it & the blue section on the right of the tube is filled with oxygen (O2) which is the Low_pressure region. In case setup we initially specified 600000pa of pressure in the high-pressure region with 101325pa in the low-pressure region.  For the initial time steps, we can observe the high-pressure region & low-pressure region have been contained in separate sections of the tube successfully by a virtual diaphragm which indicates we have used our "Events" feature correctly. Also, we can refer to the pressure plot to verify the same, from the pressure plot we observe a constant plot until 0.0011s after which it varies in a cyclic manner as the flow furthers. The varying of the plot starts as soon as the virtual diaphragm is opened i.e. after 0.001 seconds as we have mentioned in the case set up, during this time period the diaphragm will open & we observe the mixing of the two entities inside the shock tube. We observe at time step=25 the diaphragm breaks, a shock, slip surface, and expansion waves are formed and propagate through the tube. A formation of a normal shock wave, which propagates along the length of the tube until the end of the driven section which is observed at time step=50. Then it gets reflected by the wall and hence returns back, which can be noted as the mass fraction of N2 gas retracts back, and the wave furthers towards the left end. The same is depicted in the pressure contour when the right part of the tube is observed precisely. We observe when the shock wave is created, the expansion or rarefaction waves are created at the same instant. These expansion waves move in the direction opposite to that of the shock wave with the speed of sound causing continuous or smooth changes in the state of the gas. Then they get reflected by the end wall of the driver side and cause the expansion along the length of the shock tube as they move. [Note: The expansion and other changes in the state of gas happen in the direction opposite to the direction of propagation of expansion wave]. Then again, as the shock wave, after reflection from the left end, travels through the length of the tube causing the gases moving leftward to move rightward to an extent. Meanwhile, the contact surface has moved through a distance from its initial position due to the motion of the gases and the influence of pressure. Based on the above inferences the shock tube could be divided into 4 regions as shown below: 

                                                                            

                                                                                               Fig:- Shock tube shortly after diaphragm has burst.

• • Region 1: It is the region between the shock and the right end of the wall, where the gas is at rest and the shock wave has not affected it yet.
  • Region 2: It is the region between the shock wave and the contact surface, which is affected by the initial shock wave.
  • Region 3: It is the region between the contact surface and the tail of the expansion fan, which is under expansion & the gas is moving in the direction of the shock.
  • Region 4: It is the region between the head of the expansion wave & the left end of the wall, where the gas is at rest & yet to be affected by the expansion wave.  
   

C. Velocity Contour Variation Inside The Shock Tube

                                                                                                                                                  Initial Timestep

                                                      

                                                                       Time Step = 25                                                                                                                       Time Step = 50

 

                                                                       Time Step = 75                                                                                                                    Time Step = 100

 

                                                                     Time Step = 125                                                                                                                   Time Step = 150

  

                                                                    Time Step = 175                                                                                                                   Time Step = 200

  

                                                                                                                     Final Timestep

                                   

                                                                

At the initial time step for velocity Magnitude, the High_pressure region with nitrogen (N2) & the Low_pressure region with oxygen (O2) is static indicating that there is no visible movement throughout the tube crosssection as the virtual diaphragm is still closed. Once the diaphragm is opened after 0.001 seconds we start observing velocity variation throughout the shock tube cross-section. As we can observe from the animation the variation is very drastic throughout the domain. At time step=25, after the diaphragm is opened we observe the initial development of the velocity profile propagating on both sides of the tube. As the time steps increased the drastic variation in velocity is attributed to the formation of shock waves. However we observe there is no considerable change in velocity of the gases near the ends, which denotes that the gases are at rest at these locations. Throughout the entire time steps, we observe the velocity profile to be propagating only in the center of the tube.

 

D. Temperature Contour Variation Inside The Shock Tube

                                                                                                                                                 Initial Timestep

                                   

                                                                         Time Step = 25                                                                                                                   Time Step = 50

 

                                                                          Time Step = 75                                                                                                                   Time Step = 100

 

                                                                           Time Step = 125                                                                                                                   Time Step = 150

 

                                                                         Time Step = 175                                                                                                                   Time Step = 200

 

                                                                                                                                    Final Timestep

                                  

                                                                   

                                                                

                                                               

 

At the initial time step for temperature, in the case setup we initially specified 300 K of temperature for both the high-pressure & low-pressure regions. We can observe this from the animation & the temperature plot as in the initial time steps there is no variation & we observe a straight line. As the diaphragm is broken at 0.001 seconds we start observing a cyclic motion of temperature variation across the cross-section of the tube. There is a smooth variation in a particular section and a drastic variation in some sections. The contact surface and its motion can be identified clearly from the above contour by means of the temperature discontinuity prevailing in the contour in the left section of the tube. As the shock wave causes a drastic increase in the temperature as it passes by, the contact surface is found to have completely different projections of temperature when compared to the other regions. 

Conclusion                                                                                                                                                                                                                                                                                

We were able to set up a transient shock tube simulation by successfully using the "Event" feature from Converge which acted as a virtual diaphragm, we also compared the results with 2 distinct base mesh grid size & concluded that decreasing the cell size will generate more number of cells inside a geometry which allows us to capture the physics more accurately.