## 2 Methodology  
A three-dimensional model of the streamlined body was established using SOLID-WORKS program. Dimensions of the streamlined body model were as follows: radius 0.02 m, length 0.12 m (Figure 1a,b). Simulation was carried out using the CFX module from ANSYS Workbench 15, which applies the Reynolds averaged equations of motion [17]:  
\[\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial x_{j}}\big{(}\rho U_{ j}\big{)}=0 \tag{1}\]  
\[\frac{\partial\rho U_{i}}{\partial t}+\frac{\partial}{\partial x_{j}}\big{(} \rho U_{i}U_{j}\big{)}=-\frac{\partial p}{\partial x_{i}}+\frac{\partial}{ \partial x_{j}}\big{(}\tau_{ij}-\rho\overline{u_{i}u_{j}}\big{)}+S_{M} \tag{2}\]  
where \(\rho\) is the density, \(U\) is the vector of velocity, \(p\) is the pressure, \(S_{M}\) is the momentum source, \(\tau\) is the molecular stress tensor, and \(\rho\overline{u_{i}u_{j}}\) are the Reynolds stresses. The Reynolds stresses need to be modelled by additional equations of known quantities in order to achieve \"closure\". The equations used to close the system define the type of turbulence model [17]. The \(k-\varepsilon\) turbulence model was used in this study. The \(k-\varepsilon\) turbulence model [17] assumes that the turbulence viscosity is linked to the turbulence kinetic energy and dissipation via the relation:  
\[\mu_{t}=C_{\mu}\rho\frac{k^{2}}{\varepsilon} \tag{3}\]  
where \(C_{\mu}\) is a constant, \(k\) is the turbulence kinetic energy, and \(\varepsilon\) is the turbulence dissipation. The values of \(k\) and \(\varepsilon\) come directly from the differential transport equations [17] for the turbulence kinetic energy and turbulence dissipation rate:  
\[\frac{\partial(\rho\epsilon)}{\partial t}+\frac{\partial}{\partial x_{j}} \big{(}\rho U_{j}k\big{)}=\frac{\partial}{\partial x_{j}}\Bigg{[}\bigg{(}\mu +\frac{\mu_{t}}{\sigma_{k}}\bigg{)}\frac{\partial\epsilon}{\partial x_{j}} \Bigg{]}+P_{k}-\rho\varepsilon+P_{kb} \tag{4}\]  
\[\frac{\partial(\rho\epsilon)}{\partial t}+\frac{\partial}{\partial x_{j}} \big{(}\rho U_{j}\varepsilon\big{)}=\frac{\partial}{\partial x_{j}}\Bigg{[} \bigg{(}\mu+\frac{\mu_{t}}{\sigma_{\varepsilon}}\bigg{)}\frac{\partial \epsilon}{\partial x_{j}}\Bigg{]}+\frac{\varepsilon}{k}(C_{\epsilon 1}P_{k}-C_{\epsilon 2} \rho\varepsilon+C_{\epsilon 1}P_{eb}) \tag{5}\]  
where \(\mu\) is the dynamic viscosity, \(C_{\epsilon 1}\), \(C_{\epsilon 2}\), \(\sigma_{k}\), and \(\sigma_{\varepsilon}\) are constants. \(P_{kb}\) and \(P_{eb}\) represent the influence of the buoyancy forces, which are described below. \(P_{k}\) is the turbulence production due to viscous forces. The same modelling approach was adopted in the study as it is presented in [13]. A tetrahedron-type mesh was chosen for application of the finite element method (Figure 2). In order to evaluate the influence of the mesh density on the accuracy of results, a mesh sensitivity analysis was carried out. During the mesh sensitivity analysis, the mesh density was changed from 2 to 10 million elements. The changes in the values of the drag force were insignificant, as the mesh density varied from 6 million to 10 million elements. Near-surface method with inflation layers at the surface of the body (Figure 2) containing  
Figure 1: Dimensions of the streamlined body: (**a**) body; (**b**) stabilizer. 10 million elements was selected for the modeling in order to capture the boundary layer more accurately. The homogenous heat transfer models were chosen for the numerical 3D simulation [13]. The body was placed inside a channel, the outer walls of which were modeled as zero shear walls. Water velocity (\(v\)) values from the analyzed range were defined at the channel inlet while the static pressure (\(p_{s}\)) was set to be constant and equal to 0 Pa at the channel outlet. Water temperature was kept close to the boiling temperature; water pressure (\(p\)) was equal to 101,325 Pa.