# Direct Numerical Simulation of Turbulent Boundary Layer over Cubical Roughness Elements
## 2 Numerical Details  
The governing equation of incompressible flows can be expressed in non-dimensional form as follows:  
\[\frac{\partial\vec{u}_{i}}{\partial t}+\frac{\partial\vec{u}_{i}\vec{u}_{j}}{ \partial x_{j}}=-\frac{\partial\widetilde{p}}{\partial x_{i}}+\frac{1}{Re_{0} }\frac{\partial}{\partial x_{j}}\frac{\partial\vec{u}_{i}}{\partial x_{j}}+f_ {i}, \tag{1}\]  
\[\frac{\partial\vec{u}_{i}}{\partial x_{i}}=0, \tag{2}\]  
where \(x_{i}\) are the Cartesian coordinates, \(\vec{u}_{i}\) are corresponding raw velocities, and \(\widetilde{p}\) is the pressure. The velocity (\(\vec{u}\)) is decomposed into time- and ensemble-averaged (\(U\)) and fluctuating (\(u\)) components. The third term on the right-hand side of Equation (1) represents the momentum forcing (\(f_{i}\)). The governing equation is non-dimensionalized using \(U_{\infty}\) and \(\delta_{0}\), where \(\delta_{0}\) is the inlet boundary layer thickness. The Reynolds number \(Re_{0}\) is defined as \(U_{\infty}\delta_{0}/\
u\). The fractional step method [39] is adopted to discretize the governing equation by decoupling the pressure and velocity. The second-order Crank-Nicolson scheme is used to implicitly discretize the convection and viscous terms in time, and the second-order central finite difference scheme is employed to discretize all terms in space with a staggered grid. The discrete momentum forcing is explicitly determined in time based on the velocity at the previous time step to satisfy the no-slip boundary condition on the immersed boundary. The numerical algorithm for \(f_{i}\) is described in previous studies [8; 40]. A superposition of a Blasius velocity profile and isotropic free-stream turbulence is imposed on the inlet boundary. The free-stream turbulence is generated from the Orr-Sommerfeld and Squire modes in the wall-normal direction and from the Fourier modes in time and in the spanwise direction [41]. The turbulent intensity of the free-stream turbulence is set to 5% and superimposed up to 2\(\delta_{0}\), which induces the rapid decay of the free-stream turbulence in the downstream direction [42]. A convective boundary condition is used at the exit boundary, and the Neumann boundary condition is applied at the upper boundary. The no-slip boundary condition is applied at the bottom wall and to cubic roughness elements, and a periodic boundary condition is adopted in the spanwise direction. The sizes of the computational domain are 885\(\delta_{0}\), 100\(\delta_{0}\), and 37\(\delta_{0}\) in the streamwise, wall-normal, and spanwise directions, respectively. The number of grids in each direction is 4069 (\(x\)), 541 (\(y\)), and 385 (\(z\)). The grid spacing is uniform in the streamwise and spanwise directions, while the gird is stretched in the wall-normal direction using a hyperbolic tangent function: \(y(j)=L_{y}[1+\tanh(\alpha\big{\{}(j-1)/(N_{y}-1)-1\big{\}})/\tanh a],\ j=1,2, \cdots,N_{y}\). Here, \(\alpha\) is the constant 3.05, and \(L_{y}\) and \(N_{y}\) are the domain size and the number of grids in the wall-normal direction, respectively. The uniform spacing is the order of the Kolmogorov length \(\eta\) (i.e., \(\Delta x^{*}<10\eta^{*}\)) with a minimum \(\eta^{+}\) value of 2, which is sufficient to resolve all turbulence scales [43]. The superscript \(+\) represents the quantities normalized by wall units at \(Re_{\tau}=816\). The time step in wall units is 0.0241, and the averaging time is 1400\(\delta_{0}/U_{\infty}\). The parameters of the computational domain are summarized in Table 1. The procedure for the DNS of TBLs is described in previous studies [42; 44]. \begin{table}
\begin{tabular}{c c c c c c c c c c} \hline \(L_{y}\)/\(\delta_{0}\) & \(L_{y}\)/\(\delta_{0}\) & \(L_{z}\)/\(\delta_{0}\) & \(N_{x}\) & \(N_{y}\) & \(N_{z}\) & \(\Delta x^{*}\) & \(\Delta z^{*}\) & \(\Delta y_{\min}\)\({}^{*}\) & \(\Delta y_{\min}\)\({}^{*}\) \\ \hline
885 & 100 & 37 & 4609 & 541 & 385 & 11.92 & 5.97 & 0.16 & 0.95 \\ \hline \end{tabular}
\end{table} 
Table 1: Information of the computational domain.