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

Authors: Ignazio Licata; Elmo Benedetto
Journal: Journal of Computational Engineering
 (2014)

---

## 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.


---

## Body

## 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−8 sec−1, 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)ρ(g-dvdt)=∇·T,
where T is the tensor of tensions, v is the velocity, and g is the gravitational acceleration. By projecting the equation along the motion axle, in the hypothesis of uniform flow, we get
(2)ρmgsinθ+∂τxz∂z+∂τxy∂y=0,
where τxy and τxz 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]:
(3)τxy=τB+μ[(∂u/∂y)2+(∂u/∂z)2]η/2ε2+(∂u/∂y)2+(∂u/∂z)2∂u∂y,τxz=τB+μ[(∂u/∂y)2+(∂u/∂z)2]η/2ε2+(∂u/∂y)2+(∂u/∂z)2∂u∂z,
where ε is about 10−8 sec−1, 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:(4)τxy=τB+μ[(∂u/∂y)2+(∂u/∂z)2]η/2(∂u/∂y)2+(∂u/∂z)2∂u∂y,τxz=τB+μ[(∂u/∂y)2+(∂u/∂z)2]η/2(∂u/∂y)2+(∂u/∂z)2∂u∂z
as functions of ∂u/∂y and ∂u/∂z; that is, τxy=τxy(∂u/∂y,∂u/∂z) and τxz=τxz(∂u/∂y,∂u/∂z), defined in ℜ2-{(0,0)}. For simplicity we set h≝∂u/∂y and k≝∂u/∂z, getting therefore the two functions:
(5)τxy=τxy(h,k),τxz=τxz(h,k).We setk=ah in the neighborhood of the origin (0,0), so as to get a neighborhood of the origin varying the angular coefficient a. We note that this hypothesis is very binding for the stress tensor. In this case we obtain
(6)lim⁡(h,k)→(0,0)τxy(h,k)=lim⁡h→0τB+μ[h2+(ah)2]η/2h2+(ah)2h=00,
and it is an indeterminate form. Now we suppose that h and k tend to zero with the same speed and that they are in the first positive quadrant on the plane hok; this supposition is permissible since both the partial derivatives ∂u/∂y and ∂u/∂z are positive numbers; in fact the velocity is increasing along the directions y and z. Therefore in the neighborhood of the origin (0,0) we have
(7)h=k,h≥0,
and the previous limit becomes
(8)lim⁡(h,k)→(0,0)τxy(h,k)=lim⁡h→0+τB+μ[h2+h2]η/2h2+h2h=τB2,lim⁡(h,k)→(0,0)τxz(h,k)=lim⁡h→0+τB+μ[h2+h2]η/2h2+h2h=τB2.In practice we have modified the numerical code eliminating the parameterε and we have set the tensions equal to τB/2, when ∂u/∂y and ∂u/∂z are smaller than a tolerance value equal to 10-12. 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)c(y)=c0+ky.The fact that the concentration grows with the depth has been observed experimentally but it is very difficult to find a relationc(y) or c(y,z). We can consider a Taylor expansion about the point (0,0)(10)c(y,z)=c(0,0)+∂c(0,0)∂yy+∂c(0,0)∂zz+∂2c(0,0)∂y2y2+∂2c(0,0)∂z2z2+∂2c(0,0)∂y∂zyz+⋯,
or written in a dimensional form
(11)c(y,z)=c(0,0)+∂c(0,0)∂yyH+∂c(0,0)∂zzB+∂2c(0,0)∂y2(yH)2+∂2c(0,0)∂z2(zB)2+∂2c(0,0)∂y∂zyHzB+⋯,
where H and B 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
(12)c(0,0),∂c(0,0)∂y,∂c(0,0)∂z,…
using some conditions that the distribution of concentration has to satisfy. By truncating the series to the fourth order we have
(13)c(y,z)=a0+a1yH+a2zB+a3(yH)2+a4(zB)2+a5yHzB+a6(yH)3+a7(zB)3+a8yH(zB)2+a9(yH)2zB,
where
(14)a0=c(0,0),a1=∂c(0,0)∂y,a2=∂c(0,0)∂z,….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(y,(B/2)), the distribution c(y,z) is linear in accordance with (9). Therefore we have the following conditions:
(15)∂c(y,0)∂y=∂c(y,B)∂y=∂c(0,z)∂z=0,c(0,0)=cmax⁡,c(y,B2)=cmax⁡+(cmin⁡-cmax⁡)yH.These conditions allow finding the links among the coefficientsai. In fact we have(1)(16)∂c∂y=a1H+2a3yH2+a5zHB+3a6y2H3+a8z2HB2+2a9yzH2B⟹{∂c(y,0)∂y=a1H+2a3yH2+3a6y2H3=0∂c(y,B)∂y=a1H+2a3yH2+a51H+3a6y2H3∂c(y,B)∂y=+a81H+2a9yH2=0⟹(sinceitistrueforeveryy)⟹{a1=a3=a6=0a1+a5+a8=02a3+2a9=0a6=0,(I)(2)(17)∂c∂z=a2B+2a4zB2+a5yHB+3a7z2B3+2a8yzH2B+a9y2H2B⟹{∂c(0,z)∂z=a2B+2a4zB2+3a7z2B3=0⟹(sinceitistrueforeveryz)⟹{a2=0a4=0a7=0,(II)(3)(18)c(0,0)=cmax⁡⟹{a0=cmax⁡(III),(4)(19)c(y,B2)=cmax⁡+(cmin⁡-cmax⁡)yH,c(y,B2)=a0+a1yH+a22+a3(yH)2+a44+a52yH+a6(yH)3+a78+a84yH+a92(yH)2⟹(sinceitistrueforeveryy)⟹{a0+a22+a44+a78=cmax⁡a1+a52+a84=(cmin⁡-cmax⁡)a3+a92=0a6=0.(IV)By uniting the two systems, we obtain(20)a1=a2=a3=a4=a6=a7=0a0=cmax⁡,a5+a8=0,a3+a9=0,a52+a84=(cmin⁡-cmax⁡),
where we have eliminated all the linearly independent equations!By resolving the previous system we get(21)c(y,z)=cmax⁡+4(cmin⁡-cmax⁡)yzHB+4(cmax⁡-cmin⁡)yz2HB2.This relation has been modified introducing two new variablesHt and Bt 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 P1≡(Ht,0) and P2≡(0,(Bt/2)) is smaller than a value calledvelocity of threshold. Therefore the distribution of concentration is defined by
(22)c(y,z)=cmax⁡+4(cmin⁡-cmax⁡)yzHtBt+4(cmax⁡-cmin⁡)yz2HtBt2.


## 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 c(y,z), it is possible to integrate the problem with values of ε of the order of 10−8 sec−1 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.



---

*Source: 201958-2014-03-27.xml*