## 2 Numerical Modelling  
In the present work, the incompressible homogeneous density flow is considered, for which the continuity and Navier-Stokes equations in the framework of LES are given as [14]:  
\[\
abla\cdot\mathbf{U}=0 \tag{1}\]  
\[\partial_{t}\mathbf{\dot{U}}+(\mathbf{\dot{U}}\cdot\
abla)\mathbf{\dot{U}}+ \frac{1}{\rho}\
abla\bar{p}=\
abla\cdot(\boldsymbol{\tau}+\boldsymbol{\tau}^{ sgs}) \tag{2}\]  
where the bar symbol represents spatial filtering [14]. The variables in the above equations are the velocity vector \(\mathbf{U}\), the density \(\rho\), the pressure \(p\) and the viscous stress tensor \(\boldsymbol{\tau}\). The sub-grid stress tensor is given by \(\boldsymbol{\tau}^{sgs}=2v_{sgg}\mathbf{S}+\text{tr}(\boldsymbol{\tau}^{sgs} )/3\), where \(\mathbf{S}\) is the rate of the strain tensor of the resolved velocity field, \(S_{ij}=1/2(\partial\bar{U}_{i}/\partial x_{j}+\partial\bar{U}_{j}/\partial x_ {i})\), and \(v_{sgs}\) is the sub-grid viscosity. The trace of the sub-grid tensor is added to the pressure, resulting in the so-called modified pressure \(\bar{P}=\bar{p}-\text{tr}(\boldsymbol{\tau}^{sgs})/3\)[14]. Thus, the elements of the effective viscous tensor are defined as \(\tau_{ij}^{eff}=\tau_{ij}+\tau_{ij}^{ggs}=2(\
u+\upsilon_{ssg})S_{ij}\), where \(\
u\) is the kinematic viscosity, and Equation (2) can be written in the following form:  
\[\partial_{t}\bar{\mathbf{U}}+(\bar{\mathbf{U}}\cdot\
abla)\bar{\mathbf{U}}+ \frac{1}{\rho}\
abla\bar{P}=\
abla\cdot\mathbf{\tau}^{eff} \tag{3}\]  
As shown by [PERSON] et al. [15], the correct prediction of the jet dynamics and laminar-turbulent transition process, particularly at a low inlet turbulence intensity, is largely conditioned by the sub-grid viscosity model. It should not be over-dissipative, e.g., as the classical Smagorinsky model is, because an excess of an additional dissipation falsifies the transition mechanism. Therefore, in the present work, \(\upsilon_{sgs}\) is modelled by applying [PERSON]'s model [16], where \(\upsilon_{sgs}\) vanishes in laminar flows and pure shear regions. According to this model, the eddy viscosity \(\upsilon_{sgs}\) is computed as:  
\[\upsilon_{sgs}=C\sqrt{\frac{B_{\beta}}{\alpha_{ij}\alpha_{ij}}} \tag{4}\] \[\alpha_{ij}=\frac{\partial\overline{u}_{j}}{\partial x_{i}},\quad \beta_{kl}=\Delta^{2}\alpha_{mk}\alpha_{ml}\] (5) \[B_{\beta}=\beta_{11}\beta_{22}-\beta_{12}^{2}+\beta_{11}\beta_{ 33}-\beta_{13}^{2}+\beta_{22}\beta_{33}-\beta_{23}^{2} \tag{6}\]  
where the constant in (4) is taken as \(C=2.5\times 10^{-2}\)[16]. The filter width is computed as \(\Delta=(\Delta x\Delta y\Delta z)^{1/3}\), with \(\Delta x,\ \Delta y,\ \Delta z\) being the mesh spacings. The LES code employed in this study is called SAILOR. It is an academic high-order solver, which was previously verified in computations of isothermal, non-isothermal, excited and reacting jets [17]. It enables calculations in the Cartesian \((x,y,z)\) and cylindrical \((x,r,\phi)\) coordinates. Note that when the cylindrical coordinate system is used, all the spatial terms and model formulas are correspondingly transformed. A solution algorithm is based on a projection method for pressure-velocity coupling. The time integration is performed by a predictor-corrector (Adams-Bashforth/Adams-Moulton) approach. Equations (1) and (2) are discretised in space using the 6 th-order compact difference method [18]. The solution algorithm is formulated on a half-staggered mesh, where the velocity components are stored in the same computational nodes, while the pressure nodes are moved a half grid size from them. As shown in [19], the staggering of the pressure nodes is sufficient to ensure strong velocity-pressure coupling, which eliminates the well-known pressure oscillations occurring on collocated meshes. At the same time, compared to the full staggered approach, where both the velocity components and pressure values are stored in different locations, the interpolation between particular nodes is significantly reduced.