# Snow Load Shape Coefficients and Snow Prevention Method for Stepped Flat Roofs
## 3 Numerical Simulation
### Calculation Settings  
The study was carried out using ANSYS. In this work, a 2D model was adopted for simulation. The size of the numerical simulation model was the same as that of the model used in the field measurements, and the dimensions of the computational domain were 8000 mm (16 L) \(\times\) 3000 mm (15 H\({}_{2}\)), as shown in Figure 6. The length of the rear flow field was 5000 mm (10 L) to fully develop the snow flow. In the processes of snow drift, the snow concentration above 1 m from the snow surface was almost 0; so, the height of the calculation domain was considered to be 3000 mm (15 H\({}_{2}\)).  
The number of core grids accounts for more than 50% of the total grid. In order to maximize the quality of the grid around the model, an even grid with a height of 5 mm was adopted around the model. The mesh growth factor was 1.1 in the region far from the model, and the \(y^{+}\) value was controlled at about 80. The total number of meshes was about 30,000, the quality and orthogonality of the meshes were about 1, and the maximum aspect ratio of the meshes was less than 13.  
\begin{table}
\begin{tabular}{c c c} \hline \hline
**Parameter** & **Field Measurement** & **Numerical Simulation** \\ \hline Air density \(\rho\) (kg/m\({}^{3}\)) & 1.225 & 1.225 \\ Snow density \(\rho_{s}\) (kg/m\({}^{3}\)) & 120–300 & 150 \\ Snow particle diameter \(d_{s}\) (\(\upmu\)m) & 165 (Expected value) & 165 \\ Falling speed \(w_{f}\) (m/s) & 0.2–0.5 & 0.2 (\(h>\) 0.1 m), 0.35 (\(h\leq\) 0.1 m) \\ Threshold friction velocity \(u_{t}\) (m/s) & 0.2 & 0.2 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Parameters of the numerical simulation and field measurements.  
Figure 6: Computational domain and grid division.  
The measured logarithmic wind speed profile was imposed by Equation (1) at the inlet, and the turbulent kinetic energy and dissipation rate were obtained using Equations (12) and (13) [35], with \(C_{\mu}\), the empirical constant, equal to 0.09.  
\[k(z)=\frac{u_{ABL}^{*2}}{\sqrt{C_{\mu}}} \tag{12}\]  
\[\varepsilon(z)=\frac{u_{ABL}^{*3}}{\kappa z} \tag{13}\]  
The inlet snow phase was defined by the volume fraction, which refers to Equations (14) and (15) derived by [PERSON] et al. (2018) [34]:  
\[f_{sal}\ =\frac{\rho}{3.29\rho_{S}u_{*}}\bigg{(}1-\frac{u_{st}^{2}}{u_{*}^{2}} \bigg{)}, \tag{14}\]  
\[f_{susp}\ =f_{sal}\ \bigg{(}\frac{z}{h_{sal}}\bigg{)}^{-[w_{f}/(ku_{*})]}, \tag{15}\]  
where \(h_{sal}\) is the height of the saltation layer, which can be calculated using Equation (16) [18]:  
\[h_{sal}=1.6\frac{u_{*}^{2}}{2g}. \tag{16}\]  
The boundary of the ground and the model are set as the no-slip wall boundary, and the snow surface contour is updated in real time by using dynamic grid technology based on the transient method. The dynamic grid technology can consider the effect of the snow boundary on the airflow near the model and the trajectory of subsequent snow particles, to further improve the simulation accuracy [33; 34]. The principle of the dynamic grid is shown in Figure 7.  
For the values of the wall roughness height (\(K_{S}\)) and roughness constant (\(C_{S}\)), this work set the roughness constant \(C_{S}\) = 2.0 [34]. Based on Equation (17) proposed by [PERSON] et al. (2007) [36], the roughness height of the wall was obtained as \(K_{S}\) = 2 mm, which satisfies \(K_{S}\) < \(y_{p}\), and \(y_{p}\) is the center height of the minimum grid near the ground (\(y_{p}\) = 2.5 mm).  
\[K_{S}=\frac{9.793z_{0}}{C_{S}} \tag{17}\]  
The upper boundary of the domain adopts a symmetry boundary, and the outlet is set as the pressure outlet. The calculation settings are summarized in Table 2.  
Figure 7: Schematic diagram of the dynamic grid.  
To increase the credibility of our work, the methodological flow diagram is shown in Figure 8.