Numerical Modeling of an Impinging Jet Flow inside a Thermal Cut Kerf Using CFD and Schlieren Method

Abstract

The dynamics of high-pressure gas flow injected through a nozzle during a thermal cutting process has an important effect on cutting performance. In this study, an actual gas flow condition inside a cut kerf of a plasma cut sample was simulated by considering various geometric features of the cut kerf, such as kerf width difference and cutting length difference between the top and bottom surfaces. A prototype cut kerf shape was fabricated using a transparent material. A gas flow shadowgraph from inside the fabricated cut kerf was observed using the Schlieren method. In addition, image processing was performed on images obtained with the Schlieren method before and after gas injection, which were used to validate the numerical simulation models. The effect of turbulent viscosity in various turbulent models was studied using computational fluid dynamics analyses. The results confirmed that the k–ω turbulent model was the most suitable turbulent model for the geometry used in this study. The simulation results for this model were qualitatively consistent with the reported experimental measurements.

1. Introduction

Thermal cutting is one of the most important production processes and is widely employed in industrial processes such as welding, assembling, and riveting. Thermal cutting primarily means the use of energy in various forms to cut virtually any shape from iron and nonferrous materials out of sheets or large slabs. Among various thermal cutting processes, plasma arc cutting, and laser cutting are widely used in the metal producing and metalworking industries such as structural engineering, machine fabrication, energy equipment manufacturing, construction, aerospace, shipbuilding, nuclear power plant, battery, and automotive industries. Plasma arc cutting is a fabrication process that employs superheated, ionized gas funneled through a plasma torch to heat, melt and, ultimately, cut electrically conductive material into custom shapes and designs. On the other hand, the laser cutting process employs a highly concentrated beam of light that is focused on the cutting material to melt, burn, or vaporize the material at the point of focus. Nonetheless, both of these thermal cutting processes use a high-pressure assist gas to remove molten or vaporized material from the cut kerf [1].

The plasma cutting process is suitable for a wide range of metal materials, including structural steel, alloy steel, aluminum, and copper, and can cut through material thicknesses ranging between 5 mm to 150 mm for standard machines, although the range can vary depending on the plasma model [2,3,4]. Laser cutting is capable of precise control and high-speed processing, and is widely used in ultra-precision fields such as semiconductors. The application of high-power lasers is gradually expanding with the development of advanced technologies in fields such as automobiles, shipbuilding, cladding, and additive manufacturing processes. High-power lasers are widely used not only for welding but also for cutting structures [5]. In particular, the importance of laser cutting in nuclear decommissioning has been increasing recently because laser cutting has a smaller kerf width and higher cutting quality compared to plasma and flame cutting. Both plasma arc cutting, as well as laser cutting processes, can be widely used in multiple industries, and these cutting methods are the most applicable to thermal cutting processes.

During the thermal cutting process, the kinetic energy of the high-speed gas from the nozzle provides a mechanical force that ejects the melt and protects the nozzle from spatter [6]. The quality of the cutting can be evaluated from the kerf width, kerf angle, deposition of dross, bevel profile, and so on. The intensity of cutting speed, nozzle stand-off distance, and air pressure are key factors affecting the quality of the cut, among others [7]. The stand-off distance, which represents the allowable distance between the nozzle tip and workpiece, is mainly determined by the pattern and uniformity of the gas flow through the nozzle exit. Nevertheless, these parameters are still empirically determined [6].

For the gas flowing out of an axisymmetric nozzle to an ambient pressure surrounding of 1 bar, the jet at the exit becomes supersonic if the upstream total reservoir pressure exceeds 1.89 bar [8]. This is typically the case in most thermal cutting processes. In such conditions, the flow follows a Prandtl–Meyer expansion at the corner of the nozzle. These expansion waves are eventually reflected as compression waves from the constant pressure jet boundary. When the ratio of exit-to-ambient pressure is high, the compression waves of the same family intersect each other, and shockwaves start to appear as these waves coalesce. Further downstream of the flow, a shock develops in the form of a longitudinally curved surface of revolution. This shock is incident into the axis of the flow, and its slope and intensity continuously increase as the shock approaches the axis. Thus, regular reflection from the axis of flow is not possible and the appearance of Mach shock configurations appears [9]. This non-desirable aerodynamics phenomenon results in the deterioration of the dynamic characteristics of the gas flow.

Accordingly, the molten material removal rate, kerf width, surface roughness, and waviness are mainly affected by the behavior and pattern of the gas flow during the thermal cutting process. However, no reliable concept of the mechanisms of the processes inside the kerf during the cutting process has been proposed to date. Recording the processes under natural conditions is impossible because the cut walls are not transparent and the processes involve high temperatures and reflected radiation [10]. Research under natural conditions is confined to observations of particles leaving the cut channel and to inspecting the metal surface after it has been affected by the heat source [11]. Due to these restrictions, physical and mathematical modeling of the cutting processes is significantly important.

The authors in [12,13] presented open-source-based CFD solvers to simulate plasma cutting torches. Godinaud et al. [12] used a Godunov-type scheme to implement the mathematical model in OpenFOAM and validated the solver through a set of conical test cases and used the solver to simulate a three-dimensional plasma cutting torch. Park [13] included the plasma jet flow, the volume of fluid (VoF) method in identifying the gas-to-molten metal interface, and the phase change model for computing the melting process in the mathematical model to simulate plasma arc cutting. Zhou et al. [14,15] compared several turbulence models in their study to model swirling flow inside a highly constricted plasma cutting arc and implemented the model to study the effect of plasma-gas swirl flow on a plasma cutting arc.

In addition to the studies carried out in the plasma arc cutting process, several studies have investigated the gas flow pattern in the gas-assisted laser cutting process. The research works consist of both numerical simulations and experimental observations of gas flow characteristics. The authors in [16,17,18,19] investigated the interaction of gas flows within the cut kerf in order to study the dynamic characteristics of the exit jet at various stand-off distances during the laser cutting process. Darwish et al. investigated the effect of inlet stagnation pressure and nozzle geometry on the behavior of the gas flow [8]. They used a quasi 1D gas dynamics theory to calculate the exact-design operating conditions for three different supersonic nozzles. Then, they modeled the jet flow through these nozzles numerically and verified them experimentally, using Schlieren visualization. They reported that the exit jet was found to preserve its uniform distribution with parallel boundaries and low divergence under the exact-design operating condition, unlike what was observed for the other two conditions, especially for a nozzle with a small divergence angle.

The dynamic characteristics of the exit jet, from both conical and supersonic nozzles, have been comprehensively reviewed in [20,21,22]. Man et al. [20] theoretically analyzed and visualized the exit jet patterns, in the free stream, to illustrate how its dynamic characteristics are affected by the type and size of the nozzle. As a result, supersonic nozzles were found to operate more efficiently compared to conical ones for the high-pressure laser cutting process. However, all these studies were carried out in free stream.

Cho et al. [23] captured the gas flow as a function of the inlet pressure and the stand-off distance using a high-speed camera with the Schlieren method. The experiment was carried out by passing gas through two substrates placed parallel to each other to simulate the cut kerf slot. They used images before and after the gas injection to obtain the image intensity, which was applied to analyze the cutting gas flow for three different parts. They reported that the stronger the inlet pressure and the shorter the stand-off distance, the higher the image intensity value, and the higher the gas flow rate, respectively. In this study, the kerf channel was prepared without an inclined substrate to represent the cutting edge.

The authors in [17] studied and simulated the exit jet from a straight nozzle on an inclined substrate at various inclination angles. They claimed that the inclined substrate angle had a significant effect on the exit jet pattern within the cut kerf and had a negative effect on both the ability to remove molten materials and the cutting quality, due to the steep pressure gradient at a higher inclined angle. Man et al. [21] performed an experimental investigation with the shadow graph technique and demonstrated the effects of inlet stagnation pressure, nozzle tip to work-piece stand-off distance, cut kerf width, and thickness of the workpiece in relation to the behavior of the gas jet patterns inside a simulated kerf. Their findings suggested that the high-speed jet increases the ability to remove dross and consequently improves cutting quality. However, these studies did not provide a comprehensive numerical-experimental comparison of gas flow behavior inside the cut kerf slots.

The aim of the present study was to investigate the gas flow behavior inside the simulated thermal cut kerf numerically and experimentally, in a way that would be applicable to the plasma arc and laser cutting processes. A simulated cut kerf is prepared based on the real cut kerf dimension and the Schlieren method is used to visualize the gas flow behavior through the simulated kerf. The Schlieren results are then used to compare with the numerical simulation results carried out in the 3D numerical model designed based on the cut kerf geometry. The similarities in the numerical and experimental results can be used to validate the proposed model. However, the thermal state of the gas due to the heat source added during the thermal cutting process (especially during the plasma arc cutting process) and its effects on gas behavior are not covered in this study.

This paper presents an efficient numerical model to simulate the gas flow through a cut kerf that can be employed under variable conditions including nozzle exit diameter, stand-off distance, and variable pressure inlet. In this study, the numerical simulation of the gas flow was carried out using various turbulent models and the results were compared to find an optimal model which can be further used for the simulation of variable conditions as faced when thermal cutting a thick material. To the best of the authors’ knowledge, no published contributions carried out numerical modeling that can predict the gas flow behavior inside a cut kerf geometry in presence of the effects due to adjacent kerf walls. The proposed model provides a fairly good prediction compared with the experimental measurements, sufficient for a costly manufacturing process. The findings of this study can be helpful in thick metal cutting industries to increase productivity with potential cost reduction during the design process.

2. Experimental Procedure and Design

As discussed in the previous section, most of the studies carried out on nozzle flow were conducted for the free stream configuration, on an impinging substrate or through transparent plates without any substrate representing the cutting front. However, the actual flow in a thermal cutting process such as plasma arc cutting or laser cutting is not a free stream. Therefore, to simulate the actual behavior of the gas flow through the thermal cut kerf, gas flow experiments were carried out in a specially fabricated cut kerf assembly. The gas flow shadowgraph from inside the cut kerf was observed using the Schlieren method. The flow pattern obtained from the Schlerien method is then used to compare with the numerical simulation results to validate the proposed numerical model.

A conical nozzle with a straight-walled exit was used for this study. The exit diameter of the nozzle was 1.6 mm. The sample cut piece was 30 mm in thickness and the stand-off distance was 4 mm. These parameters were derived from the plasma cutting process and based on these parameters simulated cut kerf was fabricated for the Schlieren method. The nozzle and cut kerf geometries are shown in Figure 1. The preceding subsections discuss the process for fabricating the cut kerf for this study and the experimental setup for the Schlieren method.

2.1. Fabrication of a Cut Kerf Shape for the Schlieren Experiment

To visualize the dynamic behavior of a gas flow inside a cut kerf using the Schlieren method, a transparent media is required to fabricate the cut kerf with an exact design. Therefore, two transparent glass sheets were used as the walls of the kerf slot. This method was originally adopted by Man et al. [21]. The kerf shapes were measured from a sample cut piece of 30 mm thick steel plate. The offset distance in the length of the cut kerf on the top and bottom surface was used to model the inclined impinging surface of the cut front in the kerf as depicted in Figure 1a. The kerf widths on the top and bottom surfaces were measured and are shown in Figure 1b. The kerf model shown in Figure 1 was modeled from the actual plasma cut sample given in Figure 2a. Using these kerf width measurements, a 3D model with an exact cut kerf dimension was prepared. The model was used as an internal material to insert between two transparent glass plates. The glass plates and internal material assembly configuration are shown in Figure 2b. The internal material was then fabricated using a 3D printer. The assembly of the fabricated cut kerf is shown in Figure 2c. Unlike the experimental setup used for the kerf slot in this study, Man et al. used a kerf slot that had constant width at the top and bottom surface.

2.2. Experimental Setup for Schlieren Method

Gas flow through the nozzle into the fabricated cut kerf slot was observed using the Schlieren method. A Toepler’s Z-Schlieren assembly with two concave mirrors was configured as shown in Figure 3a for this experiment. This assembly is compact and avoids errors due to chromatic aberration caused by the optical lenses with the use of a condenser lens which focuses the light source into a point source that passes through a slit of 2 mm diameter. A gas supply with compressed air was connected to the nozzle via a pressure regulator. A constant pressure of 6 atm was maintained at the inlet of the nozzle for this experiment and was derived from the plasma arc cutting process. For the Schlieren setup, a 120 W high-power white LED lamp was used as the light source. The concave mirrors used in this setup were axially parabolic mirrors with a diameter of 200 mm and a focal length of 2030 mm and made of Pyrex material coated with aluminum. The reflectance of the mirror was more than 90%.

The light emitted by the high-power LED lamp is reflected by the first concave mirror and directed horizontally toward the second concave mirror. Parallel light is generated in the region between the two mirrors. In this region, the fabricated cut kerf assembly with nozzle is placed and the gas flow through the fabricated cut kerf is observed. The reflected light from the second mirror is focused and directed toward the high-speed camera. A knife edge is placed at the focal point of the second concave mirror to improve the contrast by blocking the deflected light. The experimental setup for the Schlieren method is given in Figure 3b.

The gas flow from the nozzle induces a change in gas pressure and the temperature in the test region, and causes a change in the refractive index of the region. These variations in pressure and temperature gradient can be captured with the help of the Schlieren setup. For imaging, Photron’s FASTCAM mini UX100 camera was used. No separate ND filter or band pass filter was used in the Schlieren configuration in this study. A shutter speed (exposure time) of 1/25,600 s and a frame rate of 1000 FPS were used to capture the images.

3. Theoretical Formulation of Turbulent Gas Flow Modeling

This section presents the mathematical framework developed to simulate the compressible turbulent gas flow that takes place in the gas-assisted thermal cutting process. The governing equations were solved using the commercial CFD solver ANSYS Fluent v19.0 which uses the finite volume method [24]. The flow regime is assumed to be steady, therefore steady-state simulations were performed. For pressure-velocity coupling, a coupled algorithm is employed. Turbulence is a three-dimensional phenomenon and a two-dimensional simulation approach often leads to non-physical results. Therefore, three-dimensional flow is considered for this study. The solution domain is subdivided into a finite number of contiguous control volumes and conservation equations are applied to every control volume.

The equations used for turbulent flows are obtained from those of the laminar flows using the time averaging procedure commonly known as Reynolds averaging. Thus, the flow is assumed to be governed by the compressible Reynolds-averaged Navier–Stokes (RANS) equations, and an appropriate turbulence model is applied for closure of the RANS equations to simulate all the averaged unsteadiness. The governing equations are solved in generalized coordinates and in conservative form.

Computation fluid dynamics (CFD) tools offer several turbulence models ranging from algebraic to linear and nonlinear two-equation turbulence models. For a simple viscous flow, an algebraic model does the job well because the turbulent viscosity is determined by a local function. Whereas for more complex viscous flow features such as shear layer and regions of separated flow, the two-equation turbulence model with second-order closure needs to be employed. Several studies related to the high-pressure nozzle flow [8,17,25] have employed two-equation turbulence models, namely k–ε and k–ω turbulence models, and have experimentally verified the prediction of shock wave structure and their position.

Darwish et al. [8] employed a standard k–ε turbulence model for numerical modeling of the gas-assisted laser cutting process at various stagnation pressures and successfully verified their numerical simulation result with the Schlieren method. The authors in [17] used the renormalization group (RNG) k–ε turbulence model to study the phenomena of shock wave that is induced by a supersonic impinging jet emanating from a straight nozzle onto a substrate with varying inclined angle and verified their numerical results with experimental visualization using shadowgraph imaging.

Balabel et al. [25] investigated five turbulence models, namely standard k–ε turbulence model, extended k–ε turbulence model, realizable k–ε turbulence model, shear stress transport (SST) k–ω turbulence model, and Reynolds stress model (RSM), over a wide range of nozzle pressure ratios to demonstrate their numerical accuracy in predicting the turbulent gas flow in rocket nozzle with complex nozzle wall geometry and found the SST k–ω turbulence model outperformed all other turbulence models. From their assessment, they concluded that the k–ε turbulence models perform well in most of the cases related to high-pressure gas flow; however, for the near-wall flow problems, the SST k–ω turbulence model gives better results.

In this study, the simulation model consists of gas flow through a narrow slot within the cut kerf. Thus, the flow experience narrow wall flow condition. Hence, four turbulence models, namely standard k–ε, realizable k–ε, RNG k–ε, and standard k–ω turbulence models were assessed in terms of their agreement with the experimental results. All of the simulations were carried out for the calculation domain as discussed in Section 3.2 and boundary conditions discussed in Section 3.3. These models were assessed based on the previous research carried out in the field of high-pressure nozzle flow and the nature of the flow in the current research.

3.1. Governing Equations

Based on the above assumptions, the governing equations to be solved in this numerical simulation include mass, momentum, energy, and turbulent equations and are given as follows:

Conservation of mass:

$$\nabla \bullet \left(\rho u\right)=0$$

(1)

Conservation of momentum:

$$\nabla \bullet \left(\rho uu\right)=-\nabla p+\nabla \bullet \left[\mu \left(\nabla u+{\left(\nabla u\right)}^T\right)\right]+\nabla \bullet \left[\lambda \left(\nabla \bullet u\right)I\right]$$

(2)

where the terms $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left(\nabla u+{\left(\nabla u\right)}^T\right)$ and ($\nabla \bullet u$) are the rate-of-strain tensor and the rate of expansion of the flow, respectively.

Conservation of energy:

$$\nabla \bullet \left[\rho u\left(e+\frac{1}{2}{u}^2\right)\right]=\nabla \bullet \left( \kappa \nabla T\right)+\nabla \bullet \left[-pu+\tau \bullet u\right]$$

(3)

Since in this study, we determined the working fluid as an ideal gas, the compressibility effect must follow the equation as follows:

$$p=\rho RT$$

(4)

3.2. Turbulence Models

As mentioned in the previous section, the turbulence models employed for the study consist of standard k–ε, realizable k–ε, RNG k–ε turbulence model, and the standard k–ω turbulence models and are listed in

Table 1. Simulation cases with different turbulent models.

Cases	Case 1	Case 2	Case 3	Case 4
Turbulence models	Realizable k–ε	Standard k–ε	RNG k–ε	Standard k–ω

below.

3.2.1. k–ε Turbulence Model

The turbulence kinetic energy (k) and the rate of dissipation of the turbulent kinetic energy (ε) are obtained from the following equations:

The turbulent kinetic energy equation:

$$\nabla \bullet \left(\rho uk\right)=\nabla \bullet \left(\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\nabla k\right)+G_k-\rho \epsilon$$

(5)

The turbulent kinetic energy dissipation rate equation:

$$\nabla \bullet \left(\rho u\epsilon \right)=\nabla \bullet \left(\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\nabla \epsilon \right)+C_{1\epsilon }\frac{\epsilon }{k}{G}_k-C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}$$

(6)

where ${\mu }_t$ is given by:

$${\mu }_t=\rho C_{\mu }\frac{k^2}{\epsilon }$$

(7)

In the standard k–ε turbulence model proposed by Launder and Spalding [26], all the model coefficients C_μ, C_1ε, and C_2ε are considered to be constant; however, the standard model is incapable of capturing the subtler relationships between the turbulent energy production and the turbulent stresses caused by the anisotropy of the normal stress, and therefore results in poor performance in cases of near-wall flow [27]. In the other two turbulence models (realizable and RNG) from the k–ε turbulence model family, these shortcomings of the standard k–ε model are solved by introducing a wall damping function for each of the coefficients C_μ, C_1ε, and C_2ε. This allows it to perform better in the near wall boundary condition, where the viscous sublayer is persistent. Detailed information about realizable and RNG k–ε models can be found in the research works of Shih et al. [28] and Yakhot et al. [29].

3.2.2. k–ω Turbulence Model

This turbulence model introduces a specific dissipation rate of kinetic energy (ω) instead of a dissipation rate of kinetic energy (ε). The realizable and RNG k–ε turbulence models perform fairly well in a variety of nozzle flow problems; however, they require empirical damping functions in the viscous sub-layer which are not accurate in the presence of an adverse pressure gradient [25]. In such complex flows where the separation and reattachment of flows occur, the k–ω model provides better results. Therefore, in this study, the k–ω turbulence model along with the above discussed three k–ε turbulence model was tested and compared.

In the k–ε turbulence model, the turbulent viscosity is calculated using Equation (7). The model coefficients C_μ, C_1ε, and C_2ε are considered to be constant and are determined empirically, whereas the k–ω turbulence model uses the empirical coefficients different from those of the k–ε turbulent model to calculate turbulent viscosity. The turbulent kinetic energy for this model is the same as shown in Equation (5), whereas the specific dissipation rate of kinetic energy (ω) is given by Equation (8). The turbulent viscosity equation for the k–ω turbulence model is given by Equation (9). Detailed information about the k–ω turbulent model can be found in the research work of Wilcox [24,30]. The coefficients and their values used in the four turbulence models are given in

Table 2. Coefficients for the adopted turbulence models.

Coefficients	Standard k–ε	RNG k–ε	Realizable k–ε	k–ω
$C_{\mu }$	1.44	0.0845	Calculated	-
$C_{1\epsilon }$	1.92	1.42	Calculated	-
$C_{2\epsilon }$	1.92	1.68	1.9	-
${\sigma }_k$	1	0.7194	1	2
${\sigma }_{\epsilon }$	1.3	0.7194	1.2	-
${\alpha }_{\infty }^{*}$	-	-	-	1
${\alpha }_{\infty }$	-	-	-	0.52
${\beta }_{\infty }^{*}$	-	-	-	0.072
${\beta }_i$	-	-	-	0.09
${\sigma }_{\omega }$	-	-	-	2

.

The turbulent kinetic energy dissipation rate equation:

$$\nabla \bullet \left(\rho u\omega \right)=\nabla \bullet \left(\left(\mu +\frac{{\mu }_t}{{\sigma }_{\omega }}\right)\nabla \omega \right)+\alpha \frac{\omega }{k}{G}_k-\rho \beta f_{\beta }{\omega }^2$$

(8)

where,

$${\mu }_t={\alpha }^{*}\frac{\rho k}{\omega }$$

(9)

3.3. Computational Domain and Boundary Condition

Figure 4 shows the computational domain which consists of a nozzle and a buffer region. The buffer region consists of a cut kerf channel and a free stream region extending 50 mm below the bottom surface of the cutting plate. Figure 4c shows the tetrahedral mesh grid of the computational domain. An element size of 0.30 mm was chosen for the inner walls of the nozzle, nozzle tip, and the cut kerf walls. The buffer zone under the lower surface plate was meshed to an element size of 1.00 mm and the openings (outlets) of the CFD model were meshed to 6.00 mm. The total number of mesh elements present was 866,070.

The initial and boundary conditions for velocity (u), pressure (p), temperature (T), and turbulence variables (k and ε) over the inlet, outlet, and walls were defined. All the boundaries at the openings of the computational domain were set as outlets. Figure 5 shows all the boundary layers as colored patches. The inlet and outlets are modeled as pressure inlet and pressure outlet, respectively. The inlet and outlet temperature were fixed to 25 °C and gauge inlet and outlet pressure were fixed to 6 atm and 0 atm, respectively. These values were chosen to reflect the experimental conditions. In compressible flows, isentropic relations for an ideal gas are applied to indicate total pressure, static pressure, and velocity at a pressure inlet boundary. No slip boundary condition was implemented in the walls.

The governing equation was solved using ANSYS FLUENT v19.0. The computations utilized a pressure-based iterative coupled algorithm for discretizing the convective transport terms. This algorithm solves the momentum and pressure-based continuity equations together. The full implicit coupling is achieved through an implicit discretization of pressure gradient terms in the momentum equations, and an implicit discretization of the face mass flux, including pressure dissipation terms. A compressible form of the Navier–Stokes equation along with the previously discussed turbulent models discretized by the second order upwind for the momentum, energy, and turbulence equations were used to simulate the phenomenon of flow pattern along the cut kerf slot and around the top and bottom surfaces of the cutting plate. The number of iterations was selected to be 10,000, but in these computations, all the models converged in under 7000 iterations. The convergence criteria were chosen as 10⁻⁴ for continuity, velocity, k, ε, and ω; and 10⁻⁶ for energy equations.

3.4. Mesh Independence

For a CFD model, it is necessary to gain an insight into the sensitivity of the model with respect to various changes to the model parameter values used. For this purpose, a mesh independence test was carried out on the developed model. The element sizes of the walls inside the cut kerf and inner wall of the nozzles were altered to the value range of 0.50 mm to 0.20 mm. The element size of the edges of the kerf on the top surface was changed from 0.20 mm to 0.05 mm to obtain a mesh model with various levels of refinement. Each refinement level was maintained at a ratio of 1.6. Altogether, four mesh models with cell numbers of 304,558, 501,755, 866,070, and 1,393,309 were prepared to carry out the mesh independence test using the k–ω turbulent model. Maximum velocity along the flow direction was used as a parameter to evaluate the dependency of the simulation results on mesh size. Figure 6 gives the dependency of maximum velocity with respect to the cell numbers in the mesh model. When using a finer grid, the basic flow structure changes little, but it increases the computational cost greatly. Therefore, the grid chosen represents a compromise between the accuracy and computational time.

Mesh Adaptation

The region at the stand-off distance exhibited maximum velocity and minimum static pressure. The mesh quality at this region determines the capability of the simulation model to accurately predict the shock wave structure. Using the solution obtained for the aforementioned mesh model, the mesh at the stand-off distance was refined using mesh adaptation. Using the pressure gradient value at the stand-off distance, the mesh elements were extracted. The mesh refinement was carried out for the extracted mesh elements. The refined mesh model after mesh adaption operation yielded additional 150,000 mesh elements near the nozzle exit. The mesh model obtained after mesh adaption was used to resolve the shock wave structure at the stand-off distance.

4. Results and Discussion

In this section, the results from the numerical simulations are presented first, followed by the Schlieren experiment results. The results are then compared and validated using an image processing technique.

4.1. Numerical Simulation Results

The velocity, Mach number distribution, and pressure fields along the symmetry plane of the computational domain were mapped as individual turbulence models to evaluate their validity with respect to the experimental results. The gas flow direction along the kerf slot was inclined at an angle with respect to the axis of the nozzle. Therefore, the distribution curves are partially shown in two sections: at the nozzle section where the gas flow direction is along the nozzle axis (Z1); and at the inclined section in the kerf slot and under the cutting plate (Z2). The two sections are shown with dotted lines in Figure 7. Line Z1 passes through the axis of the nozzle from the inlet to the top surface of the cutting plate and line Z2 starts from the top surface of the cutting plate at 0.5 mm offset from the cut kerf end and runs parallel to the inclined cut kerf surface.

Figure 8, Figure 9 and Figure 10 show the velocity, Mach number, and pressure distribution respectively for all the turbulence models tested in this study. In Figure 8a, the velocity field maps for the four turbulence models are shown. Figure 8b,c shows the velocity distribution curve of the turbulence models along the lines Z1 and Z2, respectively. From the velocity field figure, it can be seen that the flow exhibits an under-expansion characteristic, which is expected for the convergent subsonic nozzles operating at high exit pressure.

Figure 9a shows the Mach number distribution map of the four turbulence models. Mach number distribution curves along lines Z1 and Z2 for each model are given in Figure 9b and Figure 9c, respectively.

Similarly, Figure 10a shows the pressure field maps of the four turbulence models, and the pressure distribution curves along lines Z1 and Z2 for each model are given in Figure 10b and Figure 10c, respectively.

Table 3. Simulation cases with different turbulent models.

Cases (Turbulent Model)	Case 1 k–ε (Reliable)	Case 2 k–ε (Standard)	Case 3 k–ε (RNG)	Case 4 k–ω
No. of iterations	5418	3876	6848	5284
Computation time (s)	8474	2887	5412	4149

shows the number of iterations and the computational time taken for each case to converge.

The flow distribution along the nozzle and kerf slot for the tested turbulent models show uniformity and have good flow distribution. The velocity and pressure distribution curve in Figure 8b and Figure 10b shows high pressure and low-velocity flow at the inlet and as the geometry of the nozzle converges gradually, the flow becomes low pressure and high-velocity flow. This phenomenon is explained by Bernoulli’s principle. As the flow exits the nozzle, the flow velocity continues to rise until the Mach shock disk is observed. This phenomenon is reflected by the Mach number distribution curve in Figure 9b.

The pressure distribution curve shown in Figure 10c along the flow direction inside the cut kerf shows fluctuations and follows a harmonic wave pattern, which dampens as the flow propagates further downstream. The high- and low-pressure regions are formed along the flow direction and are shown with the crest and trough in the pressure curve. These fluctuations indicate the presence of a Mach shock disk along the flow propagation, as observed in the Schlieren experiment. Similarly, the velocity curves in Figure 8c exhibit a similar fluctuation pattern in the flow. It is evident from the velocity and pressure distribution curves that the fluctuations in these fields occur at corresponding positions and the number of fluctuations corresponds with each other. Figure 9c shows the presence of the high Mach number region in these low-pressure and high-velocity regions.

The flow behavior shown by the pressure and velocity curves corresponds with the simulation results for all turbulent models studied in this work. However, on close comparison of these models, the standard k–ε turbulent model has low fluctuation in the flow field, and the intensity of the flow field dampens very shortly after the gas exits from the nozzle. The shock wave intensity predicted by this model is comparatively lower than all other turbulence models. The standard k–ε model is incapable of capturing the subtler relationships between the turbulent energy production and the turbulent stresses caused by the anisotropy of the normal stress, and therefore it results in poor performance in cases of near-wall flow [27]. While they were suitable for modeling a free stream jet flow, the model used for the numerical simulation in this study consisted of a cut kerf that had an inclination angle. The flow from the nozzle was directed toward the inclined surface inside the kerf slot, which produces an adverse pressure gradient, which results in the substantial deterioration of the solution.

The RNG model, however, predicted that the shock Mach disks will occur at closer distances, compared to all other models, and the fluctuation in both the pressure and velocity fields is comparatively high. This model predicts that Mach shock disks that are formed in the low-pressure zone will be of higher intensity than those predicted by any other models. The k–ε (realizable) model and k–ω turbulent model, on the other hand, exhibited remarkable similarities in flow behavior. Both models predicted the location of the Mach shock disk at almost the same position. The extended k–ε models (RNG and realizable) were designed to overcome the shortcomings of the standard k–ε model by introducing wall damping functions for the model coefficients. This helps them to perform better in near-wall flow conditions. However, several researchers have mentioned in their work that these wall damping functions are not reliable [25,27] and can lead to wrongly predicted flow behavior. The high sensitivity of the RNG model is possibly due to the choice of a less reliable damping function for the given geometry. In addition, several studies have mentioned that the RNG turbulence model provides better performance with swirling flows, and can make wrong predictions when the flow is not rotating [31].

Unlike the RNG and realizable k–ε models, the k–ω turbulent model does not require wall dampening functions to predict the flow at the viscous sublayer near the wall, and can give a comparatively accurate prediction [25]. The results obtained from the CFD simulation were validated by comparison with the experimental results carried out using the Schlieren method, which indicates the model is comparatively stable and reliable for simulating flow inside a cut kerf, similar to this study.

From the results of the realizable k–ε model, it can be concluded that the performance of this model is similar to that of the k–ω turbulence model. However, from the computational cost point of view, the realizable k–ε turbulence model proves to be the most expensive. Table 3 shows that this model took the longest time in terms of convergence of the solution. Hence, for this research, k–ω turbulence model was selected as the most appropriate turbulent model.

4.2. Experimental Result from Schlieren Method

The gas flow visualized during the Schlieren experiment using a high-speed camera reveals the subtle unsteadiness in the Mach shock wave generation during the high-pressure nozzle flow. A stable Schlieren image was chosen among a series of images captured during the experiment and compared with the density gradient map obtained from the simulation results which is shown in Figure 11. The similarities between the Schlieren experiment result and the simulation results for the k–ε (realizable) and k–ω turbulence models can be seen during comparison. However, the results from RNG turbulent model shows higher fluctuation in the flow which was not visible in the experimental result whereas the standard k–ε turbulent model was not capable of showing all the expansions and contraction in flows due to the shock wave as seen in the Schlieren experiment result. Figure 12 shows the velocity map at the stand-off distance predicted by four turbulent models after the refinement of the model using mesh adaption. It can be seen that the standard k–ε turbulent model and the k–ε (realizable) turbulent model were not able to predict the low-velocity spot behind the Mach shock disk where usually the flow velocity can reach almost zero.

The image from the Schlieren experiment result was analyzed with an image processing technique, using the image processing toolbox in MATLAB. Figure 13 shows the image processing algorithm applied to the Schlieren image for image analysis in this study.

The outliers were removed to extract the region of interest (ROI) from the images obtained from the Schlieren experiment. Figure 14a,b shows the extracted ROI images obtained before and after gas injection during the Schlieren experiment, respectively. Gas movement data can be obtained from the differences between these images, by subtracting the former one from the latter [23]. Figure 14c shows the subtracted image. However, the output image from the subtraction operation has some noise. A morphological operation was applied to this image to suppress unwanted noise, and the gas flow image was obtained as shown in Figure 14d.

Figure 15 shows a side-by-side comparison of the images of the flow pattern predicted by the numerical simulation in the mesh model after mesh adaption using the k–ω turbulence model; the image obtained from the Schlieren experiment; and the result obtained after image processing, respectively. Similarities found in these images are marked by the dotted regions in Figure 15. In the dotted region ①, when the gas is emitted from the nozzle, it shows an under-expansion in flow. In the dotted region ②, it can be seen that a stagnation region has formed at the point where the under-expanded gas comes into contact with the top surface at the end of the cut kerf slot. The flow is then deflected along the top surface. In the dotted region ③, the Mach shock disk pattern inside the kerf can be seen. The comparison of these three regions in images from the Schlieren experiments and those obtained from numerical simulations show a good agreement in terms of flow structure.

In Figure 15a, the region where gas flow experienced a high flow speed and reached the sound barrier, Mach shock disks are observed as shown in Figure 15b. The larger the Mach number, the wider the size of the Mach disks [32]. As a result, the flow is deflected around these Mach disks, and the patterns shown in Figure 15 were observed. In the region where the velocity is seen to be on the rise, the density of air is found to be low, and as the sound barrier is reached the Mach shock disk is formed. In the region behind the Mach shock disk, the velocity becomes subsonic, and the density of the air builds up. This phenomenon can be clearly seen in Figure 15a–c. Thus, there should be a relationship between the flow velocity (Mach number) and the degree of flow expansion as seen in the image captured with the Schlieren technic.

For a closer examination, the outline of the flow along the kerf was plotted with respect to the reference line (line Z2) passing along the kerf slot. Figure 16a shows the reference line Z2 along the cut kerf and Figure 16b shows the outline of gas flow marked with a red curve using edge detection. Then, the distance from the reference line to the outline of gas flow was calculated as shown in Figure 16c.

After a series of image processing operations, the graph was synchronized with the density gradient fluctuation curve obtained from the simulation result of four tested turbulent models. In the position where the expansion and contraction were seen in the gas flow, low- and high-density regions can be observed from the density gradient map in Figure 11. This can be further validated by comparing the density gradient curve and the outline curve of the gas flow pattern observed in the Schlieren experiment. The comparison is given in Figure 17.

As can be seen in Section 1, in the regions indicated by lines ①, ③, and ⑤ where expansion in the gas flow is seen, the density gradient fluctuation curve shows a low-density region. Similarly, for regions ② and ④ where the contraction in flow pattern was seen, the density gradient curve shows a rise in the density. However, in Section 2, the plot obtained from the processed image includes a high amount of noise and could not be used for a meaningful comparison with the numerical simulation results. The resemblance in the results is consistent for all the turbulent models; however, the density gradient fluctuation shown by the standard k–ε turbulent model is not noticeably high, whereas for the RNG model, the gas flow fluctuation is comparatively higher in intensity and the positions of high- and low-density regions are predicted closer to each other as discussed in the previous section. Moreover, expansion and contraction in gas flow predicted by the RNG model are higher in number which was not visible during the Schlieren experiment. The comparison between results shows that the realizable k–ε turbulent model and k–ω turbulence model give better prediction among all the tested turbulence models.

As discussed in previous the section, the selection of the k–ω turbulence model was made based on the cost effectiveness in the performance. Therefore, further comparisons were made with the numerical simulation results predicted by the k–ω turbulence model. The pressure, Mach number, and velocity distribution curves obtained from the simulation results using the k–ω turbulence model were compared with the Schlieren experiment result as shown in Figure 18.

From the figure, it can be seen that in Section 1, the crests and troughs (peak and valley) in the velocity and Mach number distribution curve show some similarity with the curve drawn at the outline of gas flow distribution extracted from Schlieren results. However, in contrast to the velocity and Mach number distribution curve, the pressure distribution curve shows the opposite characteristics. Specifically, at positions indicated by lines ①, ③, and ⑤, the velocity curve, Mach number curve, and the curve extracted from the outline of gas flow show the crest, whereas the pressure curve shows the trough. The crest in the velocity curve is due to the Mach shock disk as observed in the Mach number distribution curve. Whereas at the positions indicated by lines ② and ④, the trough can be observed in velocity distribution, Mach number distribution, and gas flow outline curve and vice versa for the pressure curve. In conclusion, the similarity between the calculated velocity distribution from the numerical simulation and the results of the Schlieren experiment for the position of Mach disk distribution shows a high degree of similarities and can be applied to validate the proposed model.

5. Conclusions

The gas flow behavior inside the cut kerf slot affects the cutting quality and the performance of the gas-assisted thermal cutting process to a great extent. Therefore, studying the dynamic behavior of the flow inside the kerf can help operators effectively determine the optimum operating conditions. A kerf slot was measured using an actual sample which was cut using a thermal cutting process. Then, a representative kerf slot was fabricated per the measurement using transparent glass and a 3D-printed part. The fabricated cut kerf was used in Schlieren experiments to visualize the gas flow through the kerf slot.

In addition to the experimental study, a numerical study was carried out to investigate the gas flow in the kerf slot. For precise measurement of such gas flow dynamics, accurate numerical modeling is necessary. In this research, several Reynolds-averaged Navier–Stokes (RANS) turbulence models were used for numerical modeling of the gas flow pattern inside cut kerf: the realizable k–ε model, the standard k–ε model, the RNG k–ε model, and the k–ω model. Using these models, a CFD analysis was conducted, and the results were presented in this work. The numerical results revealed that the k–ω model and the realizable k–ε model gave the best results compared with the other models, for predicting shock waves position and the separation points. The superior performance of these models may be attributed to the formulation of these models. The realizable k–ε model incorporates some near-wall turbulence anisotropy with wall damping functions, while the k–ω turbulence model uses a specific dissipation rate of turbulent kinetic energy and hence can perform well in near-wall conditions without introducing a wall damping function. Nevertheless, the k–ω model showed a great advantage over the realizable k–ε model in computational cost. Further analysis after refining the mesh resolution in stand-off distance, the k–ω turbulent model outperformed the realizable k–ε model in modeling the shockwave structure accurately.

The validation of the proposed model was carried out by comparing the numerical simulation results with the results obtained from the Schlieren experiment. The results of the Schlieren experiment were consistent to a great extent with the predicted results from the proposed model. The position of shock waves and the separation points of the flow predicted by the proposed model strongly resembled the experimental results. The dampening of the Mach shock disk due to repeated contact with the inclined kerf slot was remarkably similar in the two results. An image processing technique was used to further support the validation of the proposed model. The plot representing the outline of the gas flow obtained from the Schlieren image showed good resemblance to the velocity, Mach number, and pressure distribution curves from the numerical simulation results, validating the proposed model for use in the study of gas flow inside the cut kerf slot.