Simulation of a snow avalanche model test using computational fluid dynamics

Abstarct

Here we apply a two-phase flow model to simulate snow avalanche motion. Flowing snow is modeled as a Bingham fluid. Small-scale laboratory experiments were conducted using a rotating viscometer to validate the constitutive model. Experimental scale model test-runs were also performed, and run-out distances and impact pressures obtained in the model tests were reproduced using the two-phase flow model. Comparisons revealed that the simulated results were strongly dependent on the basal friction angle and calculation mesh size. Although the method reproduced the laboratory model results quite well, constitutive aspects of the basal surface remained unsolved. Future research will need to incorporate a numerical technique to handle the basal boundary, such as a boundary-fitted coordinate technique.

Introduction

In 1918, a large snow avalanche struck the village of Mitsumata, Niigata Prefecture, Japan. The avalanche was approximately 400m wide and 300m long and killed 158 residents (Izumi, 1998). Since the 1930s, Japanese researchers have studied snow avalanches in the hope of minimizing avalanche damage, but financial constraints have often led to budget cuts for public works. Consequently, there remains a high snow avalanche risk.

Numerical simulations can be useful tools to estimate damage from snow avalanches. Some of the available models are AVAL-1D (Christen and others, 2002), VARA (Natale and others, 1994), MN2D (Naaim and others, 2002), TITAN2D (Patra and others, 2005) and RAMMS (Christen and others, 2010). These tools have been mainly used in Switzerland, Austria and Italy and can predict the run-out distance and flow velocity of a snow avalanche in two- and three-dimensional (2-D and 3-D) terrain. The computational demands are low. However, it is difficult for these tools to determine the vertical velocity distribution and pressure under conditions of a dynamic flow because they apply depth-averaged equations.

Other models are based on a different numerical framework. Lang and Martinelli (1979) and Lang and others (1979, 1980) proposed the AVALNCH program code, in which the numerical algorithm is based on the finite-difference method, and the Navier–Stokes equation is solved under 2-D conditions. Gauer and others (2006) developed the ANSYS CFX4 numerical tool, which is also based on a full 2-D Eulerian framework. This model uses the Bingham fluid model to describe the complex flow behavior of snow and has been used to reproduce laboratory experiments of submarine slides. The main advantage of these methods is their ability to predict vertical velocity distribution and pressure. However, in 3-D terrain, these methods are computationally demanding.

The same kind of numerical framework has been used to simulate slope disasters such as landslides, soil avalanches and slope failure. Moriguchi and others (2005) developed and applied a numerical method to simulate a real slope disaster. This model describes soil flow behavior using the Bingham fluid model, for which the shear strength is obtained using the Mohr–Coulomb failure criterion. Thus, this method can describe soil flow behavior based on typical soil parameters, such as the internal friction angle and cohesion. Modifications to the method have introduced a surface-capturing technique. Generally, it is difficult to achieve sharp free surface and volume conservation simultaneously in two-phase flow simulations based on Eulerian-type frameworks. To solve this shortcoming, Xiao and others (2005) developed the ‘tangent of hyperbola for interface capturing’ (THINC) method. The application of THINC enables simulation of sharp surfaces and conservation of the total volume of flow material. The modified method has been used to reproduce the impact force of dry sand with high accuracy (Moriguchi, 2009).

The objective of the current study was to validate the applicability of this method to the simulation of snow avalanches. First, a small-scale laboratory experiment was carried out using a rotating viscometer. The viscous characteristics of highly fluidized snow were investigated to ascertain whether the Bingham model can describe fluidized snow. Next, the model tests were carried out and snow avalanche run-out distances and impact forces were observed. Then the model tests were reproduced in the numerical simulations. Simulated results were compared with experimental results, and the effectiveness of the proposed method was assessed.

Laboratory Experiments

Snow conditions

Laboratory experiments were conducted using snow stocked in a room kept at –20˚C. This room was made available by the National Research Institute for Earth Science and Disaster Prevention (NIED), Nagaoka, Japan. Table 1 lists the snow classifications (JCMA, 1988). All the snow used in this experiment was classified as ‘lightly compacted snow’, and all experiments were carried out in the same low-temperature room.

Table 1. Classification of snow quality by density of snow (kgm^–3)

Viscosity of fluidized snow

The proposed numerical method assumes that the material is a Bingham fluid. Thus it was first necessary to confirm whether fluidized snow can be treated as a Bingham fluid. Figure 1 shows the schematic relationship between a Bingham fluid shear strain rate and shear stress. Here a rotating viscometer (Fig. 2) was used to investigate the viscous behavior of fluidized snow. Snow was placed in the cylindrical container and air was blown from the bottom of the container at a constant speed, following the procedure reported by Nishimura (1991). As shown in Figure 2, the torque of the rotating viscometer was generated by a set weight, so it was possible to maintain a constant torque. Under these conditions, we measured how long it took the fin to complete 100 revolutions (t ₁₀₀). The relation between torque and the inverse of t ₁₀₀ can be used as a substitute for the relation between shear stress and the shear strain rate. Experiments were conducted using nine different weights (0.69, 0.78, 0.98, 1.47, 1.66, 2.64, 2.94, 3.63, 4.87 N) and three kinds of material (vegetable oil, paint and snow). Figure 3 presents the results. Vegetable oil is known to be a Newton fluid, and paint is known to be a Bingham fluid, as also shown by our experimental results. In addition, the experiments confirmed that the viscous behavior of fluidized snow is similar to that in the Bingham model. Therefore, fluidized snow can be treated as a Bingham fluid. Casassa and others (1991) demonstrated that the shear strength of snow can be described using the Mohr–Coulomb failure criterion. Therefore, the proposed numerical method can describe the viscous behavior of fluidized snow.

Fig. 1. Schematic relationship between a Bingham fluid shear strain rate and shear stress.

Fig. 2. Image of the rotating viscometer.

Fig. 3. Relation between the torque and 1/t ₁₀₀.

Snow Avalanche Model Testing

A series of model tests was conducted to measure the run-out distance and impact force of snow avalanches. Tests involved a model slope (Fig. 4). Figure 5 presents the initial position of the snow mass. A triangle prism box was set at the upper section of the model slope; this box was filled with snow, and a snow flow was initiated by opening its door. The model slope was 800 cm long and 80 cm wide, and the bottom of the model slope was coated with frozen snow. If the generated snow produced snow blocks during an avalanche, it would be difficult to maintain reproducibility. To prevent snow blocks from occurring, the snow was first placed in the box and then dropped to the slope. Model testing included two kinds of tests, one measuring the run-out distance and the other measuring the impact force. Each of these involved four test cases, which were generated by changing the initial weight of the snow mass (441, 588, 736 and 883 N). For each snow weight, five trials were conducted to verify reproducibility. In all test cases, the thickness of the snow deposited on the model slope was measured, as was the initial volume of the snow mass. Table 2 lists the initial density values, which were almost identical in all cases. For impact force measurement testing, the measurement equipment was installed on the lower section of the model slope (Fig. 6). The surface of the equipment consisted of three boards, and three load cells (LCN-A-5kN: KYOWA) were set on the back of the boards to measure impact force. The boards were hung and pulled backward tightly using ropes to prevent backlash. Impact force was measured when the snow avalanche struck the equipment. The specifications of the load cells were as follows: rating capacity 5 kN, calibration coefficient 0.001251 ×10^–6 kN, rated output (RO) 2 mV V^–1 (4000×10^–6: ε) ±0.3%, nonlinearity ±0.15% RO and hysteresis ±0.10% RO.

Fig. 4. View of the model slope.

Fig. 5. Schematic side view of the model slope and the initial position of the snow mass.

Fig. 6. Measurement equipment of impact force.

Table 2. Range of initial density of snow mass (kgm^–3)

Results of Model Testing

Run-out distance measurement testing

Figure 7 shows the relationship between run-out distance and snow weight for all five trials. Run-out distances lengthened as snow weight increased. Figure 8 shows the observed velocities at different positions; the angles indicate the slope angle at each section. Flow velocities began to decrease after the snow reached a section where the slope angle was <30˚, and velocities decreased rapidly as snow reached the horizontal section.

Fig. 7. Relation between snow weight and the run-out distance for all five trials.

Fig. 8. Observed velocities at different positions; angles indicate slope angle at each section.

Impact force measurement testing

Figure 9 shows the time history of avalanche impact force during impact force measurement testing. Peak values of impact force did not occur when snow weight was 411 N, but did occur at heavier weights. In addition, impact force decreased instantaneously just after the peak. Figure 10 shows the relationship between snow weight and maximum impact force. Figure 11 shows the behavior of the snow avalanche as it made impact with the measuring equipment. The photo confirms that the front part of the snow avalanche jumped upward at impact; the decrease in impact force after the peak was caused by this jump. The results of the five trials did not vary, indicating that the impact force measurement had high reproducibility.

Fig. 9. Time history of avalanche impact force during impact force measurement testing.

Fig. 10. Relationship between snow weight and maximum impact force.

Fig. 11. Behavior of snow avalanche as it made impact with the measuring equipment.

Numerical Method

Snow modeling

Fluidized snow was assumed to be a Bingham fluid, which can be described using the relationship between shear stress and the shear strain rate. By introducing the Mohr–Coulomb failure criterion as the yield shear strength, the following equation is obtained:

(1)

where τ is the shear stress, η ₀ is the viscosity coefficient after yield, is the shear strain rate, τ _y is the yield shear strength, c is the cohesion, ϕ is the internal friction angle and p is the hydrodynamic pressure. The numerical simulation can consider an equivalent Newtonian viscosity; in case of a Newtonian fluid, the viscosity coefficient can be obtained by dividing the shear stress by the shear strain rate. The viscosity coefficient should be constant. Similarly, the equivalent Newtonian viscosity of a Bingham fluid, η ^', can

be obtained by dividing Equation (1) by the shear strain rate:

(2)

Equation (2) shows that the value of η ^' depends on the pressure and the shear strain rate. This indicates that the value of η ^' changes in both time and space. Therefore, the value of η ^' is updated at each time-step and in each calculation cell. The equation also shows that the equivalent viscosity becomes infinite as the shear strain reduces to zero. To avoid this singularity, a maximum value of the equivalent viscosity, η _max, is defined:

(3)

We used very large values for η _max (e.g. 10^˚ Pa s) because the value should theoretically be infinite. The shear strain rate shown in Equation (3) is calculated as the secondary invariant of the shear strain-rate tensor.

Governing equations

The snow is assumed to be an incompressible fluid. The following are used as governing equations:

(4)

(5)

where ui is the velocity vector, ρ is the total mass density of the snow and gi is the gravity acceleration vector. Equation (4) is the linear momentum conservation law. Equation (5) is the equation of continuity. In case of a Newtonian fluid, the viscosity coefficient is constant and its spatial derivative is zero. However, the equivalent viscosity, η ^', has a spatial derivative, so Equation (4) can be used to incorporate the spatial derivative of η ^'.

Non-advection terms, such as the pressure term, viscous term and gravity term, can be discretized using the finite-difference method. The advection term can be solved using a confined interpolation profile (CIP; Yabe and Aoki, 1991). In CIP, interpolation functions (e.g. cubic function) are used to generate a profile between calculation grids. The interpolation function moves with the velocity at each calculation grid to obtain the value at the next time-step. The key function of CIP is its ability to derive the interpolation function: the advection equations of a value and its spatial derivative are solved simultaneously. This enables two constrained conditions at each calculation grid, meaning that it is possible to build a third-order interpolation function using information at only two calculation grids.

The Poisson equation for pressure is solved implicitly, and this implicit procedure can be applied to solve the viscous term. As discussed above, equivalent viscosity, η ^', depends on the shear strain rate and calculations must be able to handle very large values. Therefore, it is necessary to use an implicit time integration scheme for the viscous term (Moriguchi and others, 2005).

The following equation can be used to capture the free surface of snow:

(6)

where φ is the volume-of-fluid (VOF) function, initially proposed by Hirt and Nichols (1981). The VOF function is defined at each calculation grid and can have a value from 0.0 to 1.0. The value indicates the occupancy of fluid at each grid. By solving Equation (6), it is possible to define the location of the surface implicitly at each time-step. In this study, the THINC method was used to solve Equation (6), making it possible to conserve the total weight of the VOF function exactly. In addition, the THINC method can maintain the shape of the fluid interface even after many calculation time-steps. These advantages are quite important for two-phase flow simulations.

Effect of basal friction angle

The proposed numerical method can also describe the effect of basal friction (Moriguchi and others, 2010). Generally, the non-slip boundary condition is described by setting the velocity vector in the opposite direction in a virtual calculation domain, as shown in Figure 12a. In contrast, Figure 12b shows that the slip boundary condition can be

Fig. 12. Boundary conditions of velocity at bottom surface: (a) non-slip boundary condition; (b) slip boundary condition.

expressed by setting the velocity vector at the same intensity at the virtual calculation point. To accommodate the effects of basal friction, parameter α, a reduction coefficient of basal friction, is applied. The velocity at the virtual calculation point, U _v, can be calculated using

(7)

where U is the velocity at a neighboring point in a calculation domain. The non-slip boundary condition and the slip boundary condition are expressed by = 1.0 and 0.0, respectively. By changing the value of from 0.0 to 1.0, it is possible to describe an arbitrary basal friction angle. In addition, can be given using

(8)

where ϕ _b is the basal friction angle. Thus, it is possible to use an arbitrary basal friction angle in the simulation.

Simulated Snow Avalanche Model Testing

Conditions of numerical analysis

Two-dimensional numerical simulations were conducted using the proposed numerical method. Two types of numerical model were prepared to simulate run-out distance measurement testing and impact force measurement testing (Fig. 13). The calculation domain was inclined to minimize calculation costs, and the incline of the model slope was specified in the form of a horizontal component of gravity. These numerical simulations only reproduced conditions of a snow weight of 883 N. Snow density was based on the results of the model testing and was set at 500 kg m^–3. Airflow was also incorporated: air density and the viscosity coefficient were set at 1.25 kg m^–3 and 2.0×10^–5 Pa s, respectively.

Fig. 13. Numerical models: (a) simulation of run-out distance measurement testing; (b) simulation of impact force measurement testing.

Table 3 shows the parameters used in this simulation. The material parameters of snow (internal friction angle ϕ and cohesion c) were based on previous research (Kamiishi and others, 2009). In the experiment, the bottom of the model slope was coated with frozen snow. Thus, the basal friction angle was considered to be lower than the internal friction angle ϕ. Two basal friction angles (5^˚ and 10˚) were used to investigate the effects of basal friction. A uniform Cartesian mesh was applied in the simulation, using a 2 cm mesh size (Δx =Δy = 2 cm), followed by 4 and 8 cm mesh sizes to assess the effects of the mesh size. Thus, the simulation involved a total of four cases, as shown in Table 4. Impact force was calculated by integrating the force acting on the wall (Fig. 13b). The impact force at each calculation grid was obtained based on hydrodynamic pressure. The Courant–Friedrichs–Lewy (CFL) condition was set at 5.0×10^–2.

Table 3. Parameters of snow and air

Table 4. Parameters used in simulation

Numerical results

Figures 14–16 present the results of the simulated run-out distance measurement tests. Figure 14 shows the simulated surface configurations of the avalanche at various times (case 2), and Figure 15 shows the simulated time history of the run-out distance. Here run-out distance is defined as the position of the front edge of the avalanche based on video images. Figure 16 shows final thickness distributions obtained during testing and simulations. Cases 1 and 2 in Figure 15 reveal that the final run-out distance lengthened as the basal friction angle decreased. Cases 2–4 in Figure 15 also show the effects of mesh size; the results in cases 3 and 4 were completely different from the experimental result. This finding indicates that simulated results are strongly dependent on mesh size, with smaller mesh sizes producing better results. Because the curved zone of the slope was described as a stepped surface due to the uniform Cartesian grid, the stepped surface had a stronger effect when a larger mesh size was used. The simulated and experimental results differed after the snow avalanche reached the curved zone. Moreover, Figure 16 reveals that the thickness distributions obtained in cases 3 and 4 did not agree with the experimental results. These results also indicate that the resolution was insufficient.

Fig. 14. Simulated surface configurations of the avalanche at various times. (Simulations of the run-out distance measurement testing, case 2.)

Fig. 15. Simulated time history of the run-out distance.

Fig. 16. Final thickness distributions obtained during testing and simulations. Results of simulation cases 1 (a), 2 (b), 3 (c) and 4 (d).

Figures 17 and 18 present the results of the simulated impact force measurement testing, with Figure 17 showing simulated surface configurations at various times (case 1). Only cases 1 and 2 were conducted during simulated impact force measurement testing. Figure 18 shows the time histories for impact force and flow velocity. Simulated impact forces had much higher peaks than did experimental impact forces. Impact forces also appeared slightly early in simulations compared to the model testing, indicating that the simulations overestimated the velocity before impact. The basal friction angles used in simulations were smaller than the real angles, but the results of the run-out distance measurement testing (Fig. 15) revealed that the simulated final run-out distances were shorter than experimental data. These two conclusions appear to be contradictory but can be explained. During simulated run-out distance measurement testing, the effect of the stepped surface was larger than that of the basal friction angle even when a small mesh size was used. In contrast, during simulated impact force measurement testing, most of the bottom surface consisted of a flat surface, and the simulations correctly revealed the effects of basal friction. Overall, the most accurate results were obtained when the basal friction was large and the mesh size was small.

Fig. 17. Simulated surface configurations at various times. (Simulation of the impact force measurement testing, case 1.)

Fig. 18. Time histories for impact force and flow velocity.

Conclusions

This study conducted model testing and simulations of snow avalanches. During model testing, run-out distances and impact forces were observed. Simulations of run-out distance measurement testing included a total of four cases by changing the basal friction angle and the mesh size. The results revealed that the final run-out distance lengthened as the basal friction angle decreased. The simulated results were strongly dependent on mesh size. The stepped surface had a strong effect even when a small mesh size was used. Simulations produced larger impact forces than those obtained during model testing. The basal friction angle used in this study was smaller than the real value.

Although the numerical method was able to predict the flow behavior of snow avalanche reasonably well, treatment of the bottom surface is still an unresolved issue. If simulations apply a very small mesh size, the effect of the stepped surface might decrease. However, calculation costs increase with the quantity of numerical mesh. Overcoming this problem will require introducing a technique to express the smooth boundary mathematically, such as the boundary-fitted coordinate technique.