Navier-Stokes Equation and Computational Scheme for Non-Newtonian Debris Flow

Abstract

This paper proposes a computational approach to debris flow model. In recent years, the theoretical activity on the classical Herschel-Bulkley model (1926) has been very intense, but it was in the early 80s that the opportunity to explore the computational model has enabled considerable progress in identifying the subclasses of applicability of different sets of boundary conditions and their approximations. Here we investigate analytically the problem of the simulation of uniform motion for heterogeneous debris flow laterally confined taking into account mainly the geological data and methodological suggestions derived from simulation with cellular automata and grid systems, in order to propose a computational scheme able to operate a compromise between “global” predictive capacities and computing effort.

1. Introduction

The mobility of granular clusters in the upper part of the mountain basins may cause a sliding similar to fluids currents. Obviously this happens in particular conditions of slope and if we have a particular solid-fluid concentration ratio. Inside the mixture that moves, there are several resistances. Despite the strong nonstationarity of currents, it is important to study the conditions of uniform or nearly uniform motion [1]. In fact, the phenomenon is so complex as to require the study of the problem in the simplest possible terms. Moreover, the equation of motion obtained in conditions of uniformity is also used in the simulations of variable motion and the uniform motion is an asymptotic condition where the natural currents tend to it in the absence of geometry variations of the contour and of supply conditions. In the case of currents in equilibrium with the bottom it was observed that the velocity profile along the wall has a zero gradient at the bottom.

The observations in correspondence with the free surface as well as the first velocity measurements made within the mixture showed that the transverse velocity profile exhibits a maximum at the centerline of the section and minimum values in correspondence with the side walls. Concentration measurements have been obtained by some authors for the case of granular mixtures devoid of fine material, in correspondence with the side walls of the channel. These measurements have shown that the concentration increases with the depth of the current, with a roughly linear trend, reaching a maximum value equal to the concentration at the bottom. In [2] the authors study the motion of a solid-liquid mixture under conditions of uniformity and, above all, they analyze the effect of the walls and the concentration distribution within the cross section. They have implemented a numerical code useful for calculating the local speed in the case of a nonhomogeneous Herschel-Bulkley fluid flowing in uniform motion within a channel. The numerical model describes the tendency of the material to sedimentation and erosion, through an iterative mechanism, for which the calculation in cells with speed lower than a critical threshold concentration is equal to that of the material at rest, while in cells with speed higher than the threshold concentration follows an increasing trend with depth.

The application of the model therefore allows simulating the different conditions that may occur for different combinations of the significant variables (bottom slope and concentration of the mixture) both in the presence of equilibrium with the erodible bed and in the presence of an interface with a nonerodible surface. The estimates of the average speed of the current obtained with the model described in [2] and those obtained assuming approximate homogeneous mixture were compared with the experimental measurements relating to a mud flow to the free surface, in uniform motion, with Herschel-Bulkley behavior. For details you can read [2]. The Herschel-Bulkley model, as is known, characterizes an entire class of strongly non-Newtonian fluid behaviors where the relationship between shear stress and shear rate is highly variable from point to point and it is impossible to define a viscosity constant. Traditionally, analytical efforts have focused on approximations able to return to the behavior described in that type of generalized Newtonian fluids. Unfortunately, the results, although mathematically flawless, have proved useful in very few cases of mainly biological interest (blood plasma). The rise of experimental mathematics has shifted the problem. The inability to find an analytical solution has been transformed in the search for useful “approximations classes” [3]. A number of papers [4–12] have fixed a solution for the Herschel-Bulkley model that has been very valuable for geological simulations. In it appears a constant, about 10⁻⁸ sec⁻¹, on which depend the tangent tensions along the direction of flow. The physical reason is fairly straightforward: this value empirically fits a typical situation of the change in viscosity and the convective motions in the natural landslides. We wondered if it is possible to generalize this solution, eliminating the dependence on epsilon, so you can get the same useful evaluations also for situations more and more distant from those in which it developed the original model. We used a suggestion that comes from the simulation with cellular automata and “grid” methods [13, 14], where the global behavior is obtained from simple local rules that describe the interactions between elementary and homogeneous volumes of fluid. This method requires careful consideration of the portion to be considered “primary” and the simulation critically depends on these choices. Change grid is equivalent to a mathematical zoom on the phenomenon, a sort of “renormalization” procedure. Therefore we worked on the idea that an analytical generalization of the model should contain a threshold value associated with the velocity of concentration, used as a key parameter on which depends much of the physical content of phenomenon and it is easily measurable. The results seem to encourage this approach, in other respects “classic” and analytical.

2. Herschel-Bulkley Solution (with Epsilon)

In the study of debris flows the use of a numerical code is crucial, since the Navier-Stokes equation is nonlinear partial differential equation and it is often impossible to solve analytically. Therefore the goal is to find solutions computationally “complacent.” This departs substantially from the traditional mathematical physics approach. One of the key physical assumptions for the calculation code allows a numerical integration of the following conservation of momentum: 𝜌𝐠−𝑑𝐯𝑑𝑡=∇⋅𝐓, (1) where 𝐓 is the tensor of tensions, 𝐯 is the velocity, and 𝐠 is the gravitational acceleration. By projecting the equation along the motion axle, in the hypothesis of uniform flow, we get 𝜌_𝑚𝑔sin𝜃+𝜕𝜏_𝑥𝑧+𝜕𝑧𝜕𝜏_𝑥𝑦𝜕𝑦=0, (2) where 𝜏_𝑥𝑦 and 𝜏_𝑥𝑧 are the tangent tensions along the flow direction (x axle).

Since the tangent tension is not defined in the region where the gradient tensor is zero, we have used a relation proposed by Bercovier and Engleman [15] and we have approximated the model of Herschel-Bulkley through the following continuous relations [7]: 𝜏_𝑥𝑦=𝜏_𝐵+𝜇𝜕𝑢/𝜕𝑦²+(𝜕𝑢/𝜕𝑧)^2𝜂/2𝜀²+𝜕𝑢/𝜕𝑦²+(𝜕𝑢/𝜕𝑧)²𝜕𝑢,𝜏𝜕𝑦_𝑥𝑧=𝜏_𝐵+𝜇𝜕𝑢/𝜕𝑦²+(𝜕𝑢/𝜕𝑧)^2𝜂/2𝜀²+𝜕𝑢/𝜕𝑦²+(𝜕𝑢/𝜕𝑧)²𝜕𝑢,𝜕𝑧 (3) where 𝜀 is about 10⁻⁸ sec⁻¹, empirical value calibrated on geophysical data of landslides, 𝜇 is generalized viscosity coefficient, and 𝜂 is the flow behavior index.

3. The Attempt to Eliminate 𝜀

As first step, the idea is to consider the following relations: 𝜏_𝑥𝑦=𝜏_𝐵+𝜇𝜕𝑢/𝜕𝑦²+(𝜕𝑢/𝜕𝑧)^2𝜂/2𝜕𝑢/𝜕𝑦²+(𝜕𝑢/𝜕𝑧)²𝜕𝑢,𝜏𝜕𝑦_𝑥𝑧=𝜏_𝐵+𝜇𝜕𝑢/𝜕𝑦²+(𝜕𝑢/𝜕𝑧)^2𝜂/2𝜕𝑢/𝜕𝑦²+(𝜕𝑢/𝜕𝑧)²𝜕𝑢𝜕𝑧 (4) as functions of 𝜕𝑢/𝜕𝑦 and 𝜕𝑢/𝜕𝑧; that is, 𝜏_𝑥𝑦=𝜏_𝑥𝑦(𝜕𝑢/𝜕𝑦,𝜕𝑢/𝜕𝑧) and 𝜏_𝑥𝑧=𝜏_𝑥𝑧(𝜕𝑢/𝜕𝑦,𝜕𝑢/𝜕𝑧), defined in ℜ²−{(0,0)}. For simplicity we set ℎ=^def𝜕𝑢/𝜕𝑦 and 𝑘=^def𝜕𝑢/𝜕𝑧, getting therefore the two functions: 𝜏_𝑥𝑦=𝜏_𝑥𝑦(ℎ,𝑘),𝜏_𝑥𝑧=𝜏_𝑥𝑧(ℎ,𝑘). (5)

We set 𝑘=𝑎ℎ in the neighborhood of the origin (0,0), so as to get a neighborhood of the origin varying the angular coefficient 𝑎. We note that this hypothesis is very binding for the stress tensor. In this case we obtain lim_{()→()ℎ,𝑘0,0}𝜏_𝑥𝑦(ℎ,𝑘)=lim_ℎ→0𝜏_𝐵+𝜇ℎ²+(𝑎ℎ)^2𝜂/2ℎ²+(𝑎ℎ)²0ℎ=0, (6) and it is an indeterminate form. Now we suppose that ℎ and 𝑘 tend to zero with the same speed and that they are in the first positive quadrant on the plane ℎ𝑜𝑘; this supposition is permissible since both the partial derivatives 𝜕𝑢/𝜕𝑦 and 𝜕𝑢/𝜕𝑧 are positive numbers; in fact the velocity is increasing along the directions 𝑦 and 𝑧. Therefore in the neighborhood of the origin (0,0) we have ℎ=𝑘,ℎ≥0, (7) and the previous limit becomes lim_{()→()ℎ,𝑘0,0}𝜏_𝑥𝑦(ℎ,𝑘)=lim_ℎ→0+𝜏_𝐵+𝜇ℎ²+ℎ^2𝜂/2ℎ²+ℎ²𝜏ℎ=^𝐵√2,lim_{()→()ℎ,𝑘0,0}𝜏_𝑥𝑧(ℎ,𝑘)=lim_ℎ→0+𝜏_𝐵+𝜇ℎ²+ℎ^2𝜂/2ℎ²+ℎ²𝜏ℎ=^𝐵√2. (8)

In practice we have modified the numerical code eliminating the parameter 𝜀 and we have set the tensions equal to 𝜏_𝐵/√2, when 𝜕𝑢/𝜕𝑦 and 𝜕𝑢/𝜕𝑧 are smaller than a tolerance value equal to 10⁻¹². Therefore we have identified the conditions that allow making assumptions of homogeneity on the behavior of the flow.

4. A New Concentration Law

In [15] the authors have used a linear distribution of concentration 𝑐𝑦=𝑐₀+𝑘𝑦. (9)

The fact that the concentration grows with the depth has been observed experimentally but it is very difficult to find a relation 𝑐(𝑦) or 𝑐(𝑦,𝑧). We can consider a Taylor expansion about the point (0,0)𝑐()+()𝑦,𝑧=𝑐0,0𝜕𝑐0,0()𝜕𝑦𝑦+𝜕𝑐0,0𝜕𝜕𝑧𝑧+²𝑐(0,0)𝜕𝑦²𝑦²+𝜕²𝑐(0,0)𝜕𝑧²𝑧²+𝜕²𝑐(0,0)𝜕𝑦𝜕𝑧𝑦𝑧+⋯, (10) or written in a dimensional form 𝑐()+()𝑦,𝑧=𝑐0,0𝜕𝑐0,0𝑦𝜕𝑦𝐻+𝜕𝑐(0,0)𝑧𝜕𝑧𝐵+𝜕²𝑐(0,0)𝜕𝑦²𝑦𝐻²+𝜕²𝑐(0,0)𝜕𝑧²𝑧𝐵²+𝜕²𝑐(0,0)𝑦𝜕𝑦𝜕𝑧𝐻𝑧𝐵+⋯, (11) where 𝐻 and 𝐵 are the height and the width of the river bed where the flow is confined. In the previous relation we have to calculate the constants 𝑐(),()0,0𝜕𝑐0,0,()𝜕𝑦𝜕𝑐0,0𝜕𝑧,… (12) using some conditions that the distribution of concentration has to satisfy. By truncating the series to the fourth order we have 𝑐𝑦,𝑧=𝑎₀+𝑎₁𝑦𝐻+𝑎₂𝑧𝐵+𝑎₃𝑦𝐻²+𝑎₄𝑧𝐵²+𝑎₅𝑦𝐻𝑧𝐵+𝑎₆𝑦𝐻³+𝑎₇𝑧𝐵³+𝑎₈𝑦𝐻𝑧𝐵²+𝑎₉𝑦𝐻²𝑧𝐵, (13) where 𝑎₀=𝑐(0,0),𝑎₁=𝜕𝑐(0,0)𝜕𝑦,𝑎₂=𝜕𝑐(0,0)𝜕𝑧,…. (14)

Global reasonable assumption is that the concentration is constant along the edges of the river bed and equal to maximum packing concentration and, instead, in the points of coordinates (𝑦,(𝐵/2)), the distribution 𝑐(𝑦,𝑧) is linear in accordance with (9). Therefore we have the following conditions: 𝜕𝑐𝑦,0=𝜕𝑦𝜕𝑐𝑦,𝐵=()𝜕𝑦𝜕𝑐0,𝑧𝑐()𝜕𝑧=0,0,0=𝑐_max,𝑐𝐵𝑦,2=𝑐_max+𝑐_min−𝑐_max𝑦𝐻. (15)

These conditions allow finding the links among the coefficients 𝑎_𝑖. In fact we have (1)𝜕𝑐=𝑎𝜕𝑦𝐻¹+2𝑎₃𝑦𝐻²+𝑎₅𝑧𝐻𝐵+3𝑎₆𝑦²𝐻³+𝑎₈𝑧²𝐻𝐵²+2𝑎₉𝑦𝑧𝐻²𝐵⟹𝜕𝑐𝑦,0=𝑎𝜕𝑦𝐻¹+2𝑎₃𝑦𝐻²+3𝑎₆𝑦²𝐻³=0𝜕𝑐𝑦,𝐵=𝑎𝜕𝑦𝐻¹+2𝑎₃𝑦𝐻²+𝑎₅1𝐻+3𝑎₆𝑦²𝐻³+𝑎₈1𝐻+2𝑎₉𝑦𝐻²⟹⟹=0sinceitistrueforevery𝑦_(I)𝑎₁=𝑎₃=𝑎₆𝑎=0₁+𝑎₅+𝑎₈=02𝑎₃+2𝑎₉𝑎=0₆=0, (16) (2)𝜕𝑐=𝑎𝜕𝑧𝐵²+2𝑎₄𝑧𝐵²+𝑎₅𝑦𝐻𝐵+3𝑎₇𝑧²𝐵³+2𝑎₈𝑦𝑧𝐻²𝐵+𝑎₉𝑦²𝐻²𝐵⟹𝜕𝑐(0,𝑧)=𝑎𝜕𝑧𝐵²+2𝑎₄𝑧𝐵²+3𝑎₇𝑧²𝐵³⟹⟹=0sinceitistrueforevery𝑧_(II)𝑎₂𝑎=0₄𝑎=0₇=0, (17) (3)𝑐(0,0)=𝑐_max⟹^(III)𝑎₀=𝑐_max, (18) (4)𝑐𝐵𝑦,2=𝑐_max+𝑐_min−𝑐_max𝑦𝐻,𝑐𝐵𝑦,2=𝑎₀+𝑎₁𝑦𝐻+𝑎₂2+𝑎₃𝑦𝐻²+𝑎₄4+𝑎₅2𝑦𝐻+𝑎₆𝑦𝐻³+𝑎₇8+𝑎₈4𝑦𝐻+𝑎₉2𝑦𝐻²⟹⟹sinceitistrueforevery𝑦_(IV)𝑎₀+𝑎₂2+𝑎₄4+𝑎₇8=𝑐_max𝑎₁+𝑎₅2+𝑎₈4=𝑐_min−𝑐_max𝑎₃+𝑎₉2𝑎=0₆=0. (19)

By uniting the two systems, we obtain 𝑎₁=𝑎₂=𝑎₃=𝑎₄=𝑎₆=𝑎₇𝑎=0₀=𝑐_max,𝑎₅+𝑎₈𝑎=0,₃+𝑎₉𝑎=0,⁵2+𝑎₈4=𝑐_min−𝑐_max, (20) where we have eliminated all the linearly independent equations!

By resolving the previous system we get 𝑐𝑦,𝑧=𝑐_max+4𝑐_min−𝑐_max𝑦𝑧𝑐𝐻𝐵+4_max−𝑐_min𝑦𝑧²𝐻𝐵². (21)

This relation has been modified introducing two new variables 𝐻_𝑡 and 𝐵_𝑡 that represent the height and the width of the river bed that are modified in the case that the flow velocity in the points of coordinates 𝑃₁≡(𝐻_𝑡,0) and 𝑃₂≡(0,(𝐵_𝑡/2)) is smaller than a value called velocity of threshold. Therefore the distribution of concentration is defined by 𝑐𝑦,𝑧=𝑐_max+4𝑐_min−𝑐_max𝑦𝑧𝐻_𝑡𝐵_𝑡+4𝑐_max−𝑐_min𝑦𝑧²𝐻_𝑡𝐵_𝑡². (22)

5. Final Results

Initially we have used the new distribution of concentration to integrate the problem through (3) without the parameter 𝜀 following the procedure described in Section 2. Unfortunately we have noted the difficulty of the code to integrate the problem in some cases and the next step has been to resolve the problem through the Bercovier and Engleman equations but using a concentration of the form suggested by relation (22). We have observed that, for very great widths of the river bed, the numerical solution of velocity, to the center, has the tendency to overlap with the velocity numerically calculated for the heterogeneous and not laterally confined flow (we call it analytic solution). By comparing with numerical solution obtained through relation (9), we can see that the new distribution is more neighbour to the analytic solution and therefore this relation gives us more interesting solutions. Besides we underline that, using the new concentration function 𝑐(𝑦,𝑧), it is possible to integrate the problem with values of 𝜀 of the order of 10⁻⁸ sec⁻¹ and instead, for such value, the linear distribution does not allow finding solutions. Let us notice that, in the simulation of the experiments, with widths of the river bed of 0.203 m we have observed that the variation of 𝜀 engraves few on the calculation of central velocity. Therefore it is possible to use values of 𝜀 that are not too small and we can get a solution in the smaller possible time. Finally we want to remember that in the precedent numerical codes there were many problems caused by the difference among the mass flows calculated in laboratory and those numerically calculated. Thanks to the new numerical code the differences are reduced, even if we are still distant from the real mass flows.

6. Conclusions

We have proposed a numerical code in the context of uniform motion for heterogeneous debris flow laterally confined. This computational scheme reflects a global description whose value critically depends on (8), which are hypothesis of “local regularity” of viscous convective motions, and it also depends on the distribution of concentration (22), which is here the variable macroscopic key. In other words, the idea is that it is precisely the variation of the resistance to deformation, not linearly proportional to the speed of angular deformation, that characterizes the typical inverse concentration distribution. “Local” hypothesis quite “elastic” allows a simple global description, easy to test computationally.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.