## 2 Numerical Modelling
### Computational Configurations  
The computational domains are shown in Figure 1. They are a simple rectangular box \(12D\times 2D\times 15D\) and cylinder \(6D\times 15D\), where \(D\) is the jet diameter. A nozzle that is present upstream of the inlet plane in real situations is not taken into account in the computations. Instead, the inlet boundary condition is specified in terms of an instantaneous velocity profile as:  
\[u(r,t)=U_{\text{mean}}(r)+u_{\text{turb}}(r,t)+u_{\text{excit}}(r,t) \tag{7}\]The velocity fluctuations \(u_{\rm turb}(r,t)\) are computed using a digital filtering method proposed by [PERSON] et al. (2018), which guarantees the properly correlated turbulent velocity field. The mean velocity profile is described by the hyperbolic tangent function as follows:  
\[U_{\rm mean}(r)=\frac{U_{0}+U_{\rm inf}}{2}-\frac{U_{0}+U_{\rm inf}}{2}\tanh \left[\frac{1}{4}\frac{R}{\theta}\left(\frac{r}{R}-\frac{R}{r}\right)\right] \tag{8}\]  
where \(U_{0}\) is the jet centerline velocity, \(U_{\rm inf}\) is the co-flow velocity, \(R=D/2\) is the radius of the jet and \(\theta\) is the momentum thickness. The co-flow \(U_{\rm inf}\) is added to compensate for the lack of a natural entrainment through the lateral walls. As mentioned in the introduction, the occurrence of the side-jets is conditioned by the appearance of relatively strong vortical structures in an initial jet region. These structures can originate from the global instability mechanism in variable density jets ([PERSON], 1976; [PERSON], 1976; [PERSON], 1976) or can be generated by an external excitation ([PERSON], 1976; [PERSON], 1976), which is also applied in the present work. It is assumed to be an axial forcing term \(u_{\rm excit}(r,t)\) in Equation (7) that is added to the axial velocity. It is defined as:  
\[u_{\rm excit}(r,t)=A\sin\left(2\pi\mathrm{St}_{\rm excit}\frac{U_{0}}{D}t\right) \tag{9}\]  
where \(A\) and \(\mathrm{St}_{\rm excit}=fD/U_{0}\) are the amplitude and the Strouhal number of the excitation, and \(f\) is the excitation frequency. At the side boundaries, the axial velocity is equal to the co-flow velocity, while the remaining velocity components are set equal to zero. The pressure is computed by the Neumann conditions at the inlet and side boundaries. At the outlet plane, the velocity components are computed using the convective boundary condition \(\partial U_{i}/\partial t+U_{C}\partial U_{i}/\partial n\), where \(U_{C}\) is a convection velocity, which is computed at every time step as a mean velocity at the outlet plane. The pressure is assumed to be constant at the outlet plane. Figure 2 shows the computational meshes used in the Cartesian and cylindrical coordinate systems in the cross-section perpendicular to the jet axis. The black circle represents the localisation of the nozzle edge. In the case of the Cartesian coordinate system, the mesh consists of 264 nodes in the axial direction and \(192\times 192\) nodes in the lateral directions. The mesh is compacted radially towards the jet region using the tangent hyperbolic function  
Figure 1: Computational domains: Cartesian (**a**) and cylindrical (**b**). and axially towards the inlet using the exponential function. The node distribution is such that across the jet, the cell sizes in the region \(y,z\leq 0.6D\) are almost uniform and equal to \(\Delta y=\Delta z=0.028D\). In the axial direction, the minimum cell size in the direct vicinity of the inlet equals \(\Delta x=0.039D\). It is worth noting that such a mesh was used by [PERSON] [22] for the numerical analysis of excited bifurcating and blooming jets. There, a comparison with the results obtained on a denser mesh (\(384\times 320\times 320\)) and experimental data for different Reynolds numbers were performed. It was shown that the mesh consisting of \(264\times 192\times 192\) nodes that is used in the present study is sufficient to correctly capture the round jet behaviour and ensures almost grid-independent solutions.