Skip to main content
Engineering LibreTexts

6.24: The Matousek (2010), (2011) Model

  • Page ID
  • \( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \) \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)\(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\) \(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\)

    6.24.1 Introduction

    The previous model, Wilson-GIW (1979) and Doron et al. (1987), are both based on spatial volumetric concentrations. The original Wilson-GIW (1979) 2 layer model assumes all solids transport takes place in the sliding bed. The Doron et al. (1987) 2 and 3 layer models assume a stationary or sliding bed with a certain concentration and velocity distribution above the bed. So here the solids transport takes place above the bed and if the bed is sliding, also in the bed. Given a certain spatial volumetric concentration, one of the outputs is the delivered volumetric concentration. To find a certain delivered concentration (this is often the operational input), one has to iterate the spatial volumetric concentration until the required delivered concentration is found.

    Matousek (2009), (2011) and Matousek & Krupicka (2010), (2011) developed a model where the delivered concentration is an input and the spatial concentration the output. This is achieved by using an erosion equation, relating the erosion rate to the bed shear stress, the so called Meyer-Peter Muller (1948) equation (MPM).

    The first version of the model was proposed for settling-slurry flow above stationary deposit ( Matousek (2009) and Matousek (2011)). It split the discharge area to the zone associated with a pipe wall and the zone associated with the top of the bed. The challenge was to evaluate bed friction and bed transport correctly in the slurry flow with a transport layer. A semi-empirical formula was suggested for ks/d (or λ12) and its shape has developed in time as our experience progressed with the flow behavior (experiments) and with the model computational stability. The last suggested version of ks/d is in Matousek & Krupicka (2014). The solids flow rate and the delivered concentration were determined using a transport formula. The formula was derived theoretically for assumed shapes of a concentration profile (linear) and a velocity profile (power law) across the transport layer (Matousek (2011)). Indeed, the resulting formula gains the form of the MPM equation, but its coefficients differ from the original empirical MPM, they are not constants. The equations for the coefficients were semi-empirically extended (Matousek (2011)) to be applicable also for transport of combined load (not only contact load, i.e. bed load).

    In the next step, the model was extended to work not only in the stationary-bed regime but also in the sliding-bed regime (Hydrotransport HT16, 2010). So, contrary to a traditional two-layer model, our model works in both regimes. The model described in this chapter is the HT16-version (2010). In 2011, the model is further modified by generalizing a determination of the bed force Fsf using the K-approach and by adding the transport-layer stress contribution to the bed force (Matousek & Krupicka (2011)).

    Figure 6.20-1 shows the two layers and the variables used. Many of the geometrical and other equations are similar or equal to the equations discussed with Wilson-GIW (1979) and Doron et al. (1987), but they are repeated here in order to have a self-containing chapter. The model is explained in 14 steps. The inputs of the model are the particle diameter d50, the liquid density ρl, the relative submerged density Rsd, the liquid viscosity νl, the pipe diameter Dp, the roughness of the pipe wall ε, the sliding friction coefficient μsf, the line speed vls, the terminal settling velocity vt, the bed concentration Cvb and the delivered volumetric concentration Cvt.

    The model in fact assumes 3 layers. The bed layer or lower layer and the upper layer above the bed. The upper layer is divided into 2 parts. A bed associated area on top of the bed and a pipe wall associated area, close to the pipe wall.

    Starting with some continuity equations, a bed height and a bed velocity, the Shields parameter is determined based on the Meyer-Peter Muller (1948) (MPM) equation, assuming the solids delivered are the result of a sliding bed and an erosion rate above the bed. This should be equal to the solids delivered based on the delivered volumetric concentration times the flow rate. Once the Shields parameter is known, the bed shear stress and the friction velocity can be determined. Based on the friction velocity and an empirical relation for the equivalent bed roughness, the Darcy Weisbach bed friction factor and bed associated radius are determined. The pipe wall associated radius follows from the bed associated radius. Now that all required parameters are known, the hydraulic gradients of the bed associated area and of the pipe wall associated area can be determined. These should be equal, since there can only be one hydraulic gradient. If the difference is significant, the bed height has to be adjusted and the calculation has to be repeated. If the difference is negligible, the equilibrium of forces on the bed is examined. If the driving force is smaller than the maximum resisting force, there is a stationary bed with a bed velocity zero. If the driving force is significantly larger than the resisting force, the bed velocity has to be increased and the calculation has to be repeated until there is an equilibrium. Finally the spatial volumetric concentration and the slip ratio can be determined.

    6.24.2 Analytical solution of the solids transport formula for shear-layer flow

    It is assumed that particles are transported exclusively within the shear layer. There is a linear concentration distribution across the shear layer. There is a power law velocity distribution in the shear layer. The velocity at the top of the shear layer is a constant factor times the friction velocity. There is a linear relationship between the bed shear stress and the thickness of the shear layer. In this book the shear layer is often referred to as sheet flow.

    A general equation for the solids flow rate in the shear layer with thickness H is:

    \[\ \mathrm{q_{s}=\int_{0}^{H} C_{v s}(y) \cdot u_{s}(y) \cdot d y}\]

    With y the vertical coordinate in the shear layer and us the local velocity in the shear layer. For the local concentration in the shear layer a linear distribution is assumed according to:

    \[\ \mathrm{C_{v s}(y)=C_{v b} \cdot \frac{H-y}{H}}\]

    So at the bottom of the shear layer (y=0) the spatial concentration equals the bed concentration and at the top of the shear layer (y=H) the spatial concentration equals zero. For the velocity in the shear layer a power law function is assumed, giving:

    \[\ \mathrm{u}_{\mathrm{s}}(\mathrm{y})=\mathrm{u}_{\mathrm{s} \mathrm{H}} \cdot\left(\frac{\mathrm{y}}{\mathrm{H}}\right)^{\mathrm{n}}\]

    So at the bottom of the shear layer the velocity is zero and at the top of the shear layer the velocity has a maximum ush. The solids flow rate can now be determined with:

    \[\ \mathrm{q}_{\mathrm{s}}=\frac{\mathrm{C}_{\mathrm{v} \mathrm{b}} \cdot \mathrm{u}_{\mathrm{s H}}}{\mathrm{H}^{\mathrm{n}+\mathrm{1}}} \cdot \int_{\mathrm{0}}^{\mathrm{H}}(\mathrm{H}-\mathrm{y}) \cdot \mathrm{y}^{\mathrm{n}} \cdot \mathrm{d} \mathrm{y}=\frac{\mathrm{C}_{\mathrm{v} \mathrm{b}} \cdot \mathrm{u}_{\mathrm{s} \mathrm{H}}}{\mathrm{H}^{\mathrm{n}+\mathrm{1}}} \cdot \int_{\mathrm{0}}^{\mathrm{H}} \mathrm{H} \cdot \mathrm{y}^{\mathrm{n}}-\mathrm{y}^{\mathrm{n}+\mathrm{1}} \cdot \mathrm{d} \mathrm{y}\]

    Integrating this over the height of the shear layer H gives:

    \[\ \begin{array}{left} \mathrm{q}_{\mathrm{s}} &=\frac{\mathrm{C}_{\mathrm{v} \mathrm{b}} \cdot \mathrm{u}_{\mathrm{s} \mathrm{H}}}{\mathrm{H}^{\mathrm{n}+1}} \cdot \int_{\mathrm{0}}^{\mathrm{H}}\left(\mathrm{H} \cdot \mathrm{y}^{\mathrm{n}}-\mathrm{y}^{\mathrm{n}+1}\right) \cdot \mathrm{d} \mathrm{y}=\frac{\mathrm{C}_{\mathrm{v} \mathrm{b}} \cdot \mathrm{u}_{\mathrm{s} \mathrm{H}}}{\mathrm{H}^{\mathrm{n}+1}} \cdot\left(\frac{\mathrm{H} \cdot \mathrm{y}^{\mathrm{n}+1}}{\mathrm{n}+\mathrm{1}}-\frac{\mathrm{y}^{\mathrm{n}+2}}{\mathrm{n}+\mathrm{2}}\right)_{\mathrm{0}}^{\mathrm{H}} \\ &=\frac{\mathrm{C}_{\mathrm{v b}} \cdot \mathrm{u}_{\mathrm{s H}}}{\mathrm{H}^{\mathrm{n}+1}}{ \cdot \frac{\mathrm{H}^{\mathrm{n}+2}}{(\mathrm{n}+\mathrm{1}) \cdot(\mathrm{n}+\mathrm{2})}}=\frac{\mathrm{C}_{\mathrm{v b}} \cdot \mathrm{u}_{\mathrm{s H}} \cdot \mathrm{H}}{(\mathrm{n}+\mathrm{1}) \cdot(\mathrm{n}+\mathrm{2})} \end{array}\]

    Now assuming that the velocity at the top of the shear layer equals the friction velocity at the bed times a constant, gives:

    \[\ \mathrm{u}_{\mathrm{s} \mathrm{H}}=\gamma \cdot \mathrm{u}_{* \mathrm{b}}\]


    \[\ \mathrm{q_{s}=\frac{C_{v b} \cdot u_{s H} \cdot H}{(n+1) \cdot(n+2)}=\frac{C_{v b} \cdot \gamma \cdot u_{* b} \cdot H}{(n+1) \cdot(n+2)}}\]

    It should be mentioned here that this solids flow has the dimension m2/s and not m/s. So multiplied by the width of the bed this gives a volume flow of m3/s.

    The Shields parameter by definition is:

    \[\ \mathrm{\theta=\frac{\mathrm{u}_{*b}^{2}}{\mathrm{R}_{\mathrm{s} d} \cdot \mathrm{g} \cdot \mathrm{d}}}\]

    The shear stress at the top of the shear layer has to be equal to the submerged weight of the shear layer times the iinternal friction coefficient:

    \[\ \mathrm{\rho_{\mathrm{l}} \cdot u_{*b}^{2}=H \cdot \frac{C_{v b}}{2} \cdot \rho_{l} \cdot R_{s d} \cdot g \cdot \tan (\varphi)}\]

    This can be written as:

    \[\ \mathrm{H}=\frac{2 \cdot \mathrm{u}_{\mathrm{*b}}^{2}}{\mathrm{C}_{\mathrm{v} \mathrm{b}} \cdot \mathrm{R}_{\mathrm{s d}} \cdot \mathrm{g} \cdot \tan (\varphi)}=\frac{2 \cdot \mathrm{u}_{\mathrm{*b}}^{2} \cdot \mathrm{d}}{\mathrm{C}_{\mathrm{v} \mathrm{b}} \cdot \mathrm{R}_{\mathrm{s d}} \cdot \mathrm{g} \cdot \mathrm{d} \cdot \tan (\varphi)}=\frac{2 \cdot \theta \cdot \mathrm{d}}{\mathrm{C}_{\mathrm{v b}} \cdot \tan (\varphi)}\]

    For the solids flow this gives:

    \[\ \mathrm{q}_{\mathrm{s}}=\frac{\gamma \cdot \mathrm{u}_{* \mathrm{b}}}{(\mathrm{n}+\mathrm{1}) \cdot(\mathrm{n}+\mathrm{2})} \cdot \frac{\mathrm{2} \cdot \mathrm{\theta} \cdot \mathrm{d}}{\tan (\varphi)}=\frac{\mathrm{2} \cdot \gamma}{(\mathrm{n}+\mathrm{1}) \cdot(\mathrm{n}+\mathrm{2}) \cdot \mathrm{t a n}(\varphi)} \cdot \sqrt{\mathrm{R}_{\mathrm{s d}} \cdot \mathrm{g} \cdot \mathrm{d}} \cdot \mathrm{d} \cdot \mathrm{\theta}^{\mathrm{1 . 5}}\]

    The total solids flow over a bed width O12 is:

    \[\ \mathrm{Q_{s}=q_{s} \cdot O_{12}=\frac{\gamma \cdot u_{* b}}{(n+1) \cdot(n+2)} \cdot \frac{2 \cdot \theta \cdot d}{\tan (\varphi)}=\frac{2 \cdot \gamma}{(n+1) \cdot(n+2) \cdot \tan (\varphi)} \cdot \sqrt{R_{s d} \cdot g \cdot d} \cdot d \cdot \theta^{1.5} \cdot O_{12}}\]

    In a non-dimensional form this can be written as:

    \[\ \mathrm{\Phi=\frac{Q_{s}}{O_{12} \cdot d \cdot \sqrt{R_{s d} \cdot g \cdot d}}=\frac{2 \cdot \gamma}{(n+1) \cdot(n+2) \cdot \tan (\varphi)} \cdot \theta^{1.5}}\]

    This is similar to the Meyer-Peter Muller (1948) (MPM) equation. For high flow velocities giving a large Shields parameter the equations are almost identical.

    Before starting the algorithm some parameters have to be determined or defined, that are constant during the iteration process:

    The particle Reynolds number Rep:

    \[\ \mathrm{Re}_{\mathrm{p}}=\frac{\mathrm{v}_{\mathrm{t}} \cdot \mathrm{d}_{\mathrm{5 0}}}{v_{\mathrm{l}}}\]

    The internal friction coefficient based on the internal friction angle φ of the bed:

    \[\ \tan (\varphi)=0.577 \quad(30^{\circ} \text{ internal friction angle for sands}) \]

    The power of the velocity distribution in the sheet flow layer n:

    \[\ \mathrm{n=1}\]

    The coefficient of the velocity at the top of the sheet flow layer γ according to Pugh & Wilson (1999):

    \[\ \gamma=9.4\]

    The Meyer-Peter Muller (1948) (MPM) equation, based on a stationary bed:

    \[\ \frac{\mathrm{Q}_{\mathrm{s}}}{\mathrm{O}_{\mathrm{1 2}} \cdot \sqrt{\mathrm{R}_{\mathrm{s d}} \cdot \mathrm{g} \cdot \mathrm{d}_{\mathrm{5} 0}^{\mathrm{3}}}}=\mathrm{\alpha} \cdot\left(\theta-\boldsymbol{\theta}_{\mathrm{c r}}\right)^{\beta}\]

    The two coefficients α and β of the MPM equation can now be determined according to:

    \[\ \begin{array}{left} \mathrm{\alpha=\frac{2 \cdot \gamma}{(n+1) \cdot(n+2) \cdot \tan (\varphi)}+\frac{58}{\operatorname{Re}_{p}^{0.62}}}\\

    The critical Shields parameter θcr can be determined using a suitable Shields diagram or equation as is discussed in Chapter 5.

    6.24.3 The Iteration Process

    Step I: The bed velocity

    The choice of the bed velocity v2 in between zero and the line speed vls.

    One can start with for example v2=0 m/sec.

    Step II: The bed height

    The choice of the bed height yb in between zero and the pipe diameter Dp.

    One can start with for example yb=0.05·Dp.

    Step III: Geometry of the discharge area

    The bed angle β can be expressed in the bed height yb by:

    \[\ \beta=\mathrm{a} \cos \left(1-\left(2 \cdot \frac{\mathrm{y}_{\mathrm{b}}}{\mathrm{D}_{\mathrm{p}}}\right)\right)\]

    The bed height yb can be expressed in the bed angle β by:

    \[\ \mathrm{y_{b}=\frac{D_{p}}{2} \cdot(1-\cos (\beta))}\]

    The cross section of the pipe Ap is:

    \[\ \mathrm{A_{p}=\frac{\pi}{4} \cdot D_{p}^{2}}\]

    The cross section of the bed A2 is now:

    \[\ \mathrm{A}_{2}=\mathrm{A}_{\mathrm{p}} \cdot \frac{(\beta-\sin (\beta) \cdot \cos (\beta))}{\pi}\]

    The cross section of the restricted area above the bed A1 is:

    \[\ \mathrm{A}_{\mathrm{1}}=\mathrm{A}_{\mathrm{p}}-\mathrm{A}_{\mathrm{2}}\]

    The length of the contact area between the upper layer (1) and the pipe wall O1 is:

    \[\ \mathrm{O}_{\mathrm{1}}=\mathrm{D}_{\mathrm{p}} \cdot(\pi-\boldsymbol{\beta})\]

    The length of the contact area between the bed (2) and the pipe wall O2 is:

    \[\ \mathrm{O}_{2}=\mathrm{D}_{\mathrm{p}} \cdot \boldsymbol{\beta}\]

    The length of the interface bed-upper layer (12) O12 is:

    \[\ \mathrm{O}_{\mathrm{1 2}}=\mathrm{D}_{\mathrm{p}} \cdot \sin (\boldsymbol{\beta})\]

    The hydraulic radius associated with the upper layer Rh1 is:

    \[\ \mathrm{R}_{\mathrm{h} 1}=\frac{\mathrm{A}_{1}}{\mathrm{O}_{1}+\mathrm{O}_{12}}\]

    Step IV: The velocity in the upper layer

    The velocity in the upper layer above the bed v1 is:

    \[\ \mathrm{v}_{1}=\left(\frac{\mathrm{v}_{\mathrm{l} \mathrm{s}} \cdot \mathrm{A}_{\mathrm{p}}-\mathrm{v}_{2} \cdot \mathrm{A}_{2}}{\mathrm{A}_{1}}\right)\]

    Step V: The relative average velocity difference between the upper layer and the bed layer

    The relative average velocity difference between the upper layer and the bed layer v12 is:

    \[\ \mathrm{v}_{12}=\mathrm{v}_{1}-\mathrm{v}_{2}\]

    Step VI: Solids flow rates

    The total solids flow rate, the delivered solids flow rate Qs is:

    \[\ \mathrm{Q}_{\mathrm{s}}=\mathrm{v}_{\mathrm{l} \mathrm{s}} \cdot \mathrm{A}_{\mathrm{p}} \cdot \mathrm{C}_{\mathrm{v} \mathrm{t}}\]

    The flow of solids in the upper layer Qs1 is:

    \[\ \mathrm{Q}_{\mathrm{s} 1}=\mathrm{v}_{\mathrm{1}} \cdot \mathrm{A}_{\mathrm{1}} \cdot \mathrm{C}_{\mathrm{v t}, \mathrm{1}}\]

    The flow of solids in the bed Qs2 is:

    \[\ \mathrm{Q}_{\mathrm{s} 2}=\mathrm{v}_{\mathrm{2}} \cdot \mathrm{A}_{\mathrm{2}} \cdot \mathrm{C}_{\mathrm{v b}}\]

    The mass continuity equation for the particle transport is:

    \[\ \begin{array}{left}\mathrm{v}_{\mathrm{l s}} \cdot \mathrm{A}_{\mathrm{p}} \cdot \mathrm{C}_{\mathrm{v t}}=\mathrm{v}_{\mathrm{1}} \cdot \mathrm{A}_{\mathrm{1}} \cdot \mathrm{C}_{\mathrm{v} \mathrm{t}, \mathrm{1}}+\mathrm{v}_{\mathrm{2}} \cdot \mathrm{A}_{2} \cdot \mathrm{C}_{\mathrm{v b}}\\
    \mathrm{Q}_{\mathrm{s} 1}=\mathrm{Q}_{\mathrm{s}}-\mathrm{Q}_{\mathrm{s} 2}=\mathrm{v}_{\mathrm{l s}} \cdot \mathrm{A}_{\mathrm{p}} \cdot \mathrm{C}_{\mathrm{v t}}-\mathrm{v}_{\mathrm{2}} \cdot \mathrm{A}_{\mathrm{2}} \cdot \mathrm{C}_{\mathrm{v} \mathrm{b}}\end{array}\]

    Intermezzo 1:

    The Shields parameter:

    \[\ \mathrm{\theta=\frac{\tau_{12}}{\rho_{\mathrm{l}} \cdot R_{s d} \cdot g \cdot d_{50}}}\]

    The Meyer-Peter Muller (1948) (MPM) equation, based on a stationary bed:

    \[\ \mathrm{\frac{Q_{s}}{O_{12} \cdot \sqrt{R_{s d} \cdot g \cdot d_{50}^{3}}}=\alpha \cdot\left(\theta-\theta_{c r}\right)^{\beta}}\]

    Step VII: Shear stress parameters

    The Shields parameter θ can be determined by the inverse MPM equation:

    \[\ \mathrm{ \theta=\left(\frac{Q_{s} \cdot \frac{v_{12}}{v_{1}}}{\alpha \cdot O_{12} \cdot \sqrt{R_{s d} \cdot g \cdot d_{50}^{3}}}\right)^{1 / \beta}+\theta_{c r}}\]

    The correction of v12/v1 is applied because the MPM equation gives the erosion rate based on a stationary bed and here is applied for a possible sliding bed. This correction is mathematically not 100% correct, but is a good approach.

    The bed shear stress \(\ \tau_{12}\) can now be determined by:

    \[\ \mathrm{\tau_{12}=\theta \cdot \rho_{\mathrm{l}} \cdot R_{\mathrm{sd}} \cdot g \cdot d_{50}}\]

    The friction velocity u* is now:

    \[\ \mathrm{u}_{*}=\sqrt{\frac{\tau_{12}}{\rho_{\mathrm{l}}}}\]

    Figure 6.24-1: The bed and wall associated areas.

    Screen Shot 2020-07-20 at 11.17.22 PM.png

    Intermezzo 2:

    The upper layer A1 can be divided in two parts. A part associated with the pipe wall A1,1 and a part associated with the bed interface A1,12, see Figure 6.24-1, so.

    \[\ \mathrm{A}_{\mathrm{1}}=\mathrm{A}_{\mathrm{1}, \mathrm{1}}+\mathrm{A}_{\mathrm{1}, \mathrm{1} \mathrm{2}}\]

    The hydraulic radius associated with the wall Rh1,1 is now:

    \[\ \mathrm{R}_{\mathrm{h} 1, \mathrm{1}}=\frac{\mathrm{A}_{1,1}}{\mathrm{O}_{1}} \Rightarrow \mathrm{R}_{\mathrm{h} 1,1} \cdot \mathrm{O}_{1}=\mathrm{A}_{1,1}\]

    The hydraulic radius associated with the bed Rh1,12 is now:

    \[\ \mathrm{R_{h 1,12}=\frac{A_{1,12}}{O_{12}} \Rightarrow R_{h 1,12} \cdot O_{12}=A_{1,12}}\]

    This gives in terms of the hydraulic gradient im:

    \[\ \mathrm{i}_{\mathrm{m}}=\lambda_{1} \cdot \frac{\mathrm{1}}{2 \cdot \mathrm{g} \cdot\left(4 \cdot \mathrm{R}_{\mathrm{h} 1,1}\right)} \cdot \mathrm{v}_{1}^{2}=\lambda_{12} \cdot \frac{\mathrm{1}}{2 \cdot \mathrm{g} \cdot\left(4 \cdot \mathrm{R}_{\mathrm{h} 1,12}\right)} \cdot \mathrm{v}_{12}^{2}\]

    Giving also:

    \[\ \mathrm{ \frac{R_{\mathrm{h} 1,1}}{R_{\mathrm{h} 1,12}}=\frac{\lambda_{1} \cdot v_{1}^{2}}{\lambda_{12} \cdot v_{12}^{2}}}\]

    The wall Darcy Weisbach friction factor λ1 can be determined with any suitable equation, like:

    \[\ \lambda_{1}=\frac{1.325}{\left(\ln \left(\frac{\varepsilon}{3.7 \cdot 4 \cdot \mathrm{R}_{\mathrm{h} 1,1}}+\frac{5.75}{\mathrm{Re}_{1}^{0.9}}\right)\right)^{2}}\]

    The bed Darcy Weisbach friction factor λ12 can be determined by:

    \[\ \sqrt{\frac{8}{\lambda_{12}}}=2.5 \cdot \ln \left(\frac{14.8 \cdot \mathrm{R}_{\mathrm{h} 1,12}}{\mathrm{k}_{\mathrm{s}}}\right)\]

    Step VIII: The bed Darcy Weisbach friction factor

    The Darcy Weisbach friction coefficient of the bed interface λ12 can now be determined from the friction velocity and the relative velocity between the upper and lower layer:

    \[\ \lambda_{12}=\frac{8 \cdot \mathrm{u}_{*}^{2}}{\mathrm{v}_{\mathrm{12}}^{2}}\]

    Step IX: The hydraulic gradient resulting from the bed friction

    The equivalent bed roughness ks can be determined with a number of equations:

    \[\ \frac{\mathrm{k}_{\mathrm{s}}}{\mathrm{d}_{\mathrm{5} 0}}=\mathrm{1. 35} \cdot\left(\left(\frac{\mathrm{R}_{\mathrm{s d}}^{2}}{\mathrm{g} \cdot v_{\mathrm{l}}}\right)^{1 / 3} \cdot \mathrm{v}_{\mathrm{t}}\right)^{0.5} \cdot \theta^{\mathrm{1.58}}\]

    The so called bed associated radius Rh1,12 now follows from:

    \[\ \begin{array}{left} \mathrm{R_{\mathrm{h} 1,12}=\frac{k_{\mathrm{s}}}{\mathrm{B}} \cdot \mathrm{e}^{\mathrm{\kappa} \cdot \sqrt{\frac{8}{\lambda_{12}}}}}\\
    \text{With: } \kappa =0.4 \text{ and } \mathrm{B}=14.8\end{array}\]

    The hydraulic gradient resulting from the bed friction im,12 is now:

    \[\ \mathrm{i}_{\mathrm{m}, \mathrm{1 2}}=\frac{\lambda_{12}}{\mathrm{R}_{\mathrm{h} 1,12}} \cdot \frac{\mathrm{v}_{12}^{\mathrm{2}}}{\mathrm{8} \cdot \mathrm{g}}\]

    Step X: The hydraulic gradient resulting from pipe wall friction

    The wall associated radius Rh1,1 can be determined once the bed associated radius Rh1,12 is known:

    \[\ \mathrm{R}_{\mathrm{h} \mathrm{1}, \mathrm{1}}=\frac{\mathrm{A}_{\mathrm{1}}-\mathrm{R}_{\mathrm{h} 1, \mathrm{1} \mathrm{2}} \cdot \mathrm{O}_{\mathrm{1} \mathrm{2}}}{\mathrm{O}_{\mathrm{1}}}\]

    The Reynolds number associated with the pipe wall is:

    \[\ \mathrm{R}_{1}=\frac{\mathrm{v}_{1} \cdot \mathrm{4} \cdot \mathrm{R}_{\mathrm{h} 1, \mathrm{1}}}{v_{\mathrm{l}}}\]

    The hydraulic gradient associated with the pipe wall im,1 can now be determined with:

    \[\ \mathrm{i}_{\mathrm{m}, \mathrm{1}}=\frac{\lambda_{\mathrm{1}}}{\mathrm{R}_{\mathrm{h} 1,1}} \cdot \frac{\mathrm{v}_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{8} \cdot \mathrm{g}}\]

    Step XI: Check the hydraulic gradients

    Both hydraulic gradients should be equal within a certain accuracy, so:

    \[\ \left|\frac{\mathrm{i}_{\mathrm{m}, \mathrm{1}}-\mathrm{i}_{\mathrm{m}, \mathrm{1} 2}}{\mathrm{i}_{\mathrm{m}, \mathrm{1}}}\right|<\text{ accuracy }\quad\text{ or }\quad\left|\frac{\mathrm{i}_{\mathrm{m}, \mathrm{1}}-\mathrm{i}_{\mathrm{m}, \mathrm{1} 2}}{\mathrm{i}_{\mathrm{m}, \mathrm{12}}}\right|< \text{ accuracy }\]

    If this is not the case, the bed height yb has to be adjusted and the iteration has to be repeated from Step II. If the two hydraulic gradients are equal within a certain accuracy, the hydraulic gradient im equals:

    \[\ \mathrm{i}_{\mathrm{m}}=\frac{\mathrm{i}_{\mathrm{m}, \mathrm{1}}+\mathrm{i}_{\mathrm{m}, \mathrm{1 2}}}{2}\]

    Step XII: Determine the driving and resisting forces on the bed

    The driving force on the bed Fdr with a pipe length ΔL equals the pressure gradient on the bed plus the shear stress at the interface times the surface of the bed:

    \[\ \mathrm{F}_{\mathrm{d r}}=\Delta \mathrm{p} \cdot \mathrm{A}_{2}+\tau_{12} \cdot \mathrm{O}_{12} \cdot \Delta \mathrm{L}=\mathrm{i}_{\mathrm{m}} \cdot \rho_{\mathrm{l}} \cdot \mathrm{g} \cdot \Delta \mathrm{L} \cdot \mathrm{A}_{2}+\tau_{12} \cdot \mathrm{O}_{12} \cdot \mathrm{\Delta} \mathrm{L}\]

    The resisting force on the bed Fsf equals the normal force between bed and pipe surface times a sliding friction coefficient, Matousek & Krupicka (2010):

    \[\ \mathrm{F_{\mathrm{sf}}=\mu_{\mathrm{sf}} \cdot \rho_{\mathrm{l}} \cdot \mathrm{R}_{\mathrm{sd}} \cdot \mathrm{g} \cdot \mathrm{C}_{\mathrm{vb}} \cdot \mathrm{A}_{\mathrm{p}} \cdot \frac{\mathrm{2} \cdot(\sin (\beta)-\beta \cdot \cos (\beta))}{\pi} \cdot \Delta \mathrm{L}}\]

    Later they added the normal force due to the weight of the sheet flow layer, Matousek & Krupicka (2011).

    \[\ \mathrm{F_{sf}=\mu_{sf}\cdot\rho_{l}\cdot R_{sd}\cdot g \cdot C_{vb}\cdot A_p \cdot \frac{2}{\pi}\cdot} \left(\begin{array}{l}\frac{\mathrm{K}-\mathrm{1}}{\mathrm{3}} \cdot \sin ^{3}(\beta)+\sin (\beta)-\frac{\mathrm{K}+\mathrm{1}}{\mathrm{2}} \beta \cdot \cos (\beta) \\ +\frac{\mathrm{K}-\mathrm{1}}{2} \cdot \cos ^{2}(\beta) \cdot \sin (\beta)\end{array}\right) \cdot \Delta \mathrm{L}\\
    +\frac{\tau_{12} \cdot \mathrm{D}_{\mathrm{p}}}{\tan (\varphi)} \cdot\left(\frac{\mathrm{K}+\mathrm{1}}{\mathrm{2}} \cdot \beta-\frac{\mathrm{K}-\mathrm{1}}{\mathrm{2}} \cdot \cos (\beta) \cdot \sin (\beta)\right) \cdot \Delta \mathrm{L}\]

    With K=1 this gives:

    \[\ \mathrm{F_{\mathrm{sf}}=\mu_{\mathrm{sf}} \cdot \rho_{\mathrm{l}} \cdot R_{\mathrm{sd}} \cdot \mathrm{g} \cdot \mathrm{C}_{\mathrm{vb}} \cdot \mathrm{A}_{\mathrm{p}} \cdot \frac{2 \cdot(\sin (\beta)-\beta \cdot \cos (\beta))}{\pi} \cdot \Delta \mathrm{L}+\frac{\tau_{12} \cdot \mathrm{D}_{\mathrm{p}} \cdot \beta}{\tan (\varphi)} \cdot \Delta \mathrm{L}}\]

    This resisting sliding friction force is based on the hydrostatic normal stress approach of Wilson et al. (1992), which is questioned by Miedema & Ramsdell (2014). They suggest the bed submerged weight approach giving:

    \[\ \mathrm{F}_{\mathrm{s f}}=\mu_{\mathrm{s f}} \cdot \rho_{\mathrm{l}} \cdot \mathrm{R}_{\mathrm{s d}} \cdot \mathrm{g} \cdot \mathrm{C}_{\mathrm{v b}} \cdot \mathrm{A}_{2}=\mu_{\mathrm{s f}} \cdot \rho_{\mathrm{l}} \cdot \mathrm{R}_{\mathrm{s d}} \cdot \mathrm{g} \cdot \mathrm{C}_{\mathrm{v b}} \cdot \mathrm{A}_{\mathrm{p}} \cdot \frac{\left(\beta^{\prime}-\sin \left(\beta^{\prime}\right) \cdot \cos \left(\beta^{\prime}\right)\right)}{\pi}\]

    Now there is a difference in approach between the two methods. Matousek & Krupicka (2010) exclude the sheet flow layer in their calculation of the bed angle β, while Miedema & Ramsdell (2014) include the sheet flow layer. This means that the Miedema & Ramsdell (2014) equation gives the weight of all solids in the pipe, while the Matousek & Krupicka (2011) equation adds the contribution of the weight of the sheet flow layer separately, which follows from the bed shear stress and is an essential part of the Matousek & Krupicka (2010) model. The normal stress exerted by the sheet flow layer on the bed has to be corrected for the hydrostatic normal stress approach.

    This gives a modified equation according to:

    \[\ \mathrm{F_{\mathrm{sf}}=\mu_{\mathrm{sf}} \cdot\left(\rho_{\mathrm{l}} \cdot R_{\mathrm{sd}} \cdot g \cdot C_{\mathrm{vb}} \cdot A_{\mathrm{p}}+\frac{\tau_{12} \cdot O_{12}}{\tan (\varphi)} \cdot \frac{\pi}{(\beta-\sin (\beta) \cdot \cos (\beta))}\right) \cdot \frac{2 \cdot(\sin (\beta)-\beta \cdot \cos (\beta))}{\pi} \cdot \Delta \mathrm{L}}\]

    Step XIII: Check the force balance on bed

    In the case of a stationary bed, the resisting force on the bed is bigger than the driving force and the velocity of the bed is zero, a reason to start the iteration with a bed velocity v2 of zero. If the driving force however is bigger than the resisting force, the bed velocity v2 should be increased and the iteration should be repeated from step 1 until there is a force equilibrium.

    Step XIV: Spatial concentration and slip ratio

    Above the bed a sheet flow layer is assumed with an average concentration of 50% of the bed concentration. The height of this sheet flow layer follows from the shear stress and the weight of the sheet flow layer giving:

    \[\ \mathrm{H}=\frac{2 \cdot \tau_{12}}{\rho_{\mathrm{l}} \cdot \mathrm{R}_{\mathrm{sd}} \cdot \mathrm{g} \cdot \mathrm{C}_{\mathrm{v b}} \cdot \tan (\varphi)}\]

    With the definition of the Shields parameter this can be rewritten to:

    \[\ \mathrm{H}=\frac{2 \cdot \theta \cdot \mathrm{d}_{50}}{\mathrm{C}_{\mathrm{v b}} \cdot \tan (\varphi)}\]

    The cross section of the sheet flow layer reduced to the bed concentration is:

    \[\ \mathrm{A_{12, \mathrm{sf}}=\frac{H}{2} \cdot O_{12}=\frac{\theta \cdot d_{50}}{C_{\mathrm{vb}} \cdot \tan (\varphi)} \cdot O_{12}}\]

    The total cross section containing solids with the bed concentration is now, only counting the solids, not the pores:

    \[\ \mathrm{A}_{\mathrm{s}}=\left(\mathrm{A}_{2}+\mathrm{A}_{\mathrm{1} 2, \mathrm{s f}}\right) \cdot \mathrm{C}_{\mathrm{v b}}=\left(\mathrm{A}_{\mathrm{p}} \cdot \frac{(\beta-\sin (\beta) \cdot \cos (\beta))}{\pi}+\frac{\theta \cdot \mathrm{d}_{50}}{\mathrm{C}_{\mathrm{v b}} \cdot \tan (\varphi)} \cdot \mathrm{O}_{12}\right) \cdot \mathrm{C}_{\mathrm{v} \mathrm{b}}\]

    This gives a spatial volumetric concentration of:

    \[\ \mathrm{C}_{\mathrm{v s}}=\frac{\mathrm{A}_{\mathrm{s}}}{\mathrm{A}_{\mathrm{p}}}=\frac{(\beta-\sin (\beta) \cdot \cos (\beta))}{\pi} \cdot \mathrm{C}_{\mathrm{v b}}+\frac{\theta \cdot \mathrm{d}_{50}}{\mathrm{A}_{\mathrm{p}} \cdot \tan (\varphi)} \cdot \mathrm{O}_{12}\]

    The slip ratio is defined in this book as the ratio of the slip velocity to the line speed ξ=vsl/vls, with:

    \[\ \mathrm{C}_{\mathrm{v} t}=\left(1-\frac{\mathrm{v}_{\mathrm{sl}}}{\mathrm{v}_{\mathrm{ls}}}\right) \cdot \mathrm{C}_{\mathrm{v s}}\]

    The slip ratio equals:

    \[\ \xi=\frac{\mathrm{v}_{\mathrm{s} \mathrm{l}}}{\mathrm{v}_{\mathrm{ls}}}=\mathrm{1}-\frac{\mathrm{C}_{\mathrm{v} \mathrm{t}}}{\mathrm{C}_{\mathrm{v s}}}\]

    Matousek & Krupicka (2010) use another formulation for the determination of the spatial volumetric concentration:

    \[\ \mathrm{C_{v s}=\frac{C_{v b} \cdot A_{2}+C_{v t, 1} \cdot A_{1}}{A_{p}}}\]

    This formulation is correct in case of a stationary bed, assuming there is hardly slip in the upper layer, so the delivered concentration and the spatial concentration are almost equal. However with a sliding bed or with significant slip in the upper layer, this approach is not very accurate and the concentration in the upper layer does not have to be equal to the initial delivered concentration.

    6.24.4 The Concentration Distribution, Including a Sliding Bed

    Matousek & Krupicka (2014) and Matousek et al. (2014) developed the 1-D SDM, 1-Dimensional Stress Distribution based Model, in order to determine the concentration distribution of pipe flow with the possibility of a stationary or sliding bed. The algorithm of this method can be found in these publications. Considerations

    There are three mechanisms that can support particles in a flowing slurry. The first one is the interaction of particles with turbulent eddies of the carrier liquid. The second one is the interaction of particles with other particles sporadic by collisions and the third one permanently in a bed. Particles that do not contribute to the bed, can be dispersed either by diffusive action of turbulent eddies or by particle-particle collisions. The different mechanisms produce different shapes of concentration profiles. It is however very well possible that the mechanisms occur at the same time in a flow, each mechanism supporting a certain fraction of the particles.

    Recently attempts have been made to use CFD to model heterogeneous slurry flows in pipes. The advantage is more insight in the internal structure of the flow, the disadvantage is that each simulation is like an experiment, so to find general trends many simulations have to be carried out. Practice still requires engineering models. Reason to develop a one dimensional model predicting the concentration along the vertical axis of a flow cross section, assuming a constant concentration in the transverse direction.

    Typical one dimensional models are the open channel models based on turbulent diffusion of the Schmidt-Rouse type. In a general form and a stationary situation, the upwards flow of particles due to diffusion equals the downwards flow due to gravity (hindered settling);

    \[\ \begin{array}{left}\mathrm{q}_{\mathrm{s}, \mathrm{u p}}(\mathrm{z})=\mathrm{q}_{\mathrm{s}, \mathrm{d o w n}}(\mathrm{z}) \Rightarrow-\varepsilon_{\mathrm{s}} \cdot \frac{\mathrm{d} \mathrm{C}_{\mathrm{v s}}(\mathrm{z})}{\mathrm{d} \mathrm{z}}=\mathrm{C}_{\mathrm{v s}}(\mathrm{z}) \cdot \mathrm{v}_{\mathrm{t}} \cdot\left(\mathrm{1}-\mathrm{C}_{\mathrm{v s}}(\mathrm{z})\right)^{\beta}\\
    \mathrm{C}_{\mathrm{v s}}(\mathrm{z}) \cdot \mathrm{v}_{\mathrm{t}} \cdot\left(\mathrm{1}-\mathrm{C}_{\mathrm{v s}}(\mathrm{z})\right)^{\beta}+\varepsilon_{\mathrm{s}} \cdot \frac{\mathrm{d} \mathrm{C}_{\mathrm{v s}}(\mathrm{z})}{\mathrm{d} \mathrm{z}}=\mathrm{0}\end{array}\]

    When the power β=1 the Hunt (1954) equation is found. Some model modifications introduce the effect of a broad PSD of a transported solid fraction (Karabelas (1977)). A key parameter of the turbulent diffusion model is the particle diffusion coefficient εs. This is modelled in different ways in different implementations of the model. Kaushal & Tomita (2013) improved their equation, including the influence of the mean particle diameter dm to pipe diameter ratio, giving for narrow graded sands:

    \[\ \begin{array}{left}\beta_{\mathrm{sm}}=1+93.77 \cdot\left(\frac{\mathrm{d}_{\mathrm{m}}}{\mathrm{D}_{\mathrm{p}}}\right) \cdot \mathrm{e}^{1.055 \cdot \sigma_{\mathrm{g}} \cdot \frac{\mathrm{C}_{\mathrm{vs}}}{\mathrm{C}_{\mathrm{vb}}}}\\
    \text{With : }\quad \sigma_{\mathrm{g}}\text{ in range 1.15-4 }\quad\text{ and }\quad 93.77 \cdot\left(\frac{\mathrm{d}_{\mathrm{m}}}{\mathrm{D}_{\mathrm{p}}}\right) \quad\text{ in range 0.125-2.5}\end{array}\]

    They stated that their model is appropriate to different flow patterns, including stationary and sliding beds. The question arises however whether coarser grains encounter more intense turbulent support in partially stratified flow. The question also arises whether higher concentrations lead to stronger turbulent support (partially stratified flow). One could ask the question whether the turbulent diffusion equation is valid in partially stratified flow in the first place. Kaushal & Tomita (2013) use the Hunt equation and the Karabelas methodology. To implement the effect of concentration on the terminal settling velocity they use the local concentration to determine the local terminal hindered settling velocity in the solution of the Karabelas method. The Hunt equation however already contains the influence of the upwards liquid flow, because there is a downwards particle flow. Using the hindered settling in the solution of the Karabelas method sort of doubles the effect of hindered settling. At least a factor 1 should be deducted from the hindered settling power to take into account the effect of the upwards liquid velocity as already present in the Hunt equation. The Kaushal & Tomita (2013) models are not capable to simulate the measured sharp change in the concentration gradient at the top of a sliding bed.

    Gillies & Shook (1994) modified the general turbulent diffusion equation to make it more suitable for stratified flows. The presence of contact load is represented by the coefficient κ:

    \[\ -\varepsilon_{\mathrm{s}} \cdot \frac{\mathrm{d} \mathrm{C}_{\mathrm{vs}}(\mathrm{z})}{\mathrm{d} \mathrm{z}}=\sqrt{1-\kappa(\mathrm{z})} \cdot \mathrm{C}_{\mathrm{vs}}(\mathrm{z}) \cdot \mathrm{v}_{\mathrm{t}} \cdot\left(1-\mathrm{C}_{\mathrm{vs}}(\mathrm{z})\right)^{\beta}\]

    This coefficient is defined as the ratio of the contact support force to the submerged weight of the particle. If this coefficient equals 1, all particles contribute to contact load, resulting in a bed. If the coefficient equals zero., all particles contribute to suspended load and the model becomes the normal turbulent diffusion model. A value of this coefficient is not proposed however. If the coefficient is z-dependent, like in the above equation, then it could serve as a free empirical parameter taking care of a perfect match between predicted and measured local concentration values, but only if the concentration profile is already known.

    Recent analysis of Matousek & Krupicka (2014) and Matousek et al. (2014) has shown that the use of turbulent diffusion models with high values of the particle diffusivity do not seem to be appropriate for stratified flows with concentration gradients of almost zero. Stratified flows with a transport layer (shear layer or sheet flow) require an almost constant concentration gradient in the transport layer and a zero concentration gradient in the sliding or stationary bed below the transport layer. It is difficult to accept that the turbulent diffusion concept with the terminal (hindered) setting velocity of the particle and the particle turbulent diffusivity as key parameters is an appropriate modeling approach to determine concentration profiles in stratified (high concentration) flows. Hybrid modelling taking contact and suspended loads into account separately is a preferred approach. Note of the author: Suspended load may behave according to a turbulent diffusion model, a stationary or sliding bed behaves according to soil mechanics with internal and external friction, which is not included in turbulent diffusion. Summary of the Model

    The abstract of Matousek & Krupicka (2014) gives a good summary of what the model is capable of:

    A one-dimensional profile of solids concentration is modeled in a cross section of partially stratified flow of slurry in a pipe. In the flow, a certain proportion of solids is transported as contact load and occupies a sliding bed and a transport layer above the bed. The rest of solid particles is transported as suspended load within and above the transport layer. A model based on a vertical distribution of shear stress is developed to predict a concentration profile in a cross section of such flow. Besides the concentration profile, the model predicts the thickness and velocity of the bed. Furthermore, the model determines a vertical position of the top of the transport layer and a position of the hydrodynamic axis of the flow. Model predictions show a satisfactory match with new experimentally determined profiles collected in slurry flows of four different fractions of glass beads in a 100-mm pipe of their laboratory loop.

    Figure 6.24-2, Figure 6.24-3, Figure 6.24-4 and Figure 6.24-5 show experimental concentration profiles as measured by Matousek & Krupicka (2014) in a Dp=0.1 m pipe with particles of d=0.44 mm and d=0.53 mm. The Limit Deposit Velocities are about 2.5 m/s and 2.8 m/s for sand. Since here the particles are spherical glass beads, the LDV values may be a bit smaller. The velocities used in the experiments were about 3.5 m/s and 4 m/s, so above the Limit Deposit Velocities. The experiments are compared with the concentration profiles determined with the DHLLDV Framework as described in chapter 7.10. A bed concentration of 60% is assumed for all cases, based on the maximum bottom concentration measured. The bed concentration should not be changed, it is a fixed value for a certain type of particles. The bottom concentration however will change depending on the spatial concentration and the line speed to LDV ratio. Two settings were used for the coefficient in the hindered settling power. A value of 2.75 for the d=0.44 mm particles and the default value of 4 (see chapter 7.10) for the d=0.53 mm particles. The d=0.9 mm particles (not shown here) also require a factor 4. This factor appears to be constant for a certain type of particle.

    From the figures it is clear that a higher line speed gives a steeper concentration profile and a smaller bottom concentration. In general the simulated concentration profiles match the experimental values well.

    Figure 6.24-2: Experiments of Matousek & Krupicka (2014) in a Dp=0.1 m pipe with d=0.44 mm particles, Cvs=0.23.

    Screen Shot 2020-07-21 at 9.20.50 AM.png

    Figure 6.24-3: Experiments of Matousek & Krupicka (2014) in a Dp=0.1 m pipe with d=0.44 mm particles, Cvs=0.40.

    Screen Shot 2020-07-21 at 9.22.07 AM.png

    Figure 6.24-4: Experiments of Matousek & Krupicka (2014) in a Dp=0.1 m pipe with d=0.53 mm particles, Cvs=0.20.

    Screen Shot 2020-07-21 at 9.24.07 AM.png

    Figure 6.24-5: Experiments of Matousek & Krupicka (2014) in a Dp=0.1 m pipe with d=0.53 mm particles, Cvs=0.34.

    Screen Shot 2020-07-21 at 9.26.06 AM.png

    6.24.5 Conclusions & Discussion

    Matousek (2009), (2011) and Matousek & Krupicka (2010), (2011) developed a new approach for 2LM and 3 LM models. Their approach is in fact a 3LM model with a bed area, a bed associated area and a pipe wall associated area. The model is based on the delivered volumetric concentration as an input and the spatial volumetric concentration, hydraulic gradient, bed height and other parameters as output, after iterations. Wilson-GIW (1979) and Doron et al. (1987) have the spatial volumetric concentration as input and the delivered volumetric concentration as output, after iterations.

    The concept of the model is very interesting, however there are still some questions about the implementation. Central in the model are two empirical equations. The Meyer-Peter Muller (1948) equation for the erosion rate and the equation of Matousek & Krupicka (2014) for the equivalent bed roughness. It should be noted, though, that the model transport equation has a theoretical background as shown in the derivation by Matousek (2011). The formula has been validated for Cvt up to almost 0.30 and the entire model for even higher Cvt values (up to almost 0.40).

    The MPM equation is designed for erosion in open channel flow with a free surface, giving the “Law of the Wall” velocity distribution. In pipe flow there will be a completely different velocity profile in the cross section of the pipe. Further it is the question whether MPM can deal with very high erosion rates and volumetric concentrations as present in slurry flow.

    The coefficients of the MPM equation are based on integration of the concentration and velocity profile in a sheet flow layer. Based on the Shields parameter found and the width of the bed interface, the cross section of the sheet flow layer can be determined. This cross section however has no relation with the so called bed associated area, which would be expected.

    The sliding friction force is determined based on the hydrostatic normal stress approach of Wilson et al. (1992). This approach may give still reasonable values for the sliding friction force up to a bed angle of 90°, but above 90° it strongly overestimates the sliding friction force. Miedema & Ramsdell (2014) investigated this force and came up with the weight approach. A result of the hydrostatic approach may be that at high concentrations the error in the sliding friction force is compensated with a to small sliding friction coefficient.

    The sliding friction resisting force misses the weight of the sheet flow layer, corrected for the hydrostatic normal stress approach. The correct equation is given.

    The relation for the spatial volumetric concentration may be correct for the case of a stationary bed, but not for the case of a sliding bed. The correct equation is given by the authors.

    The wall associated area does not contain any solids in the model, this could be an interesting extension of the model.

    There is no check to see if the sheet flow layer either touches the top of the pipe, resulting in a different concentration profile and velocity profile, or occupies the whole bed, resulting in heterogeneous flow. At very high concentrations, resulting in a high value of the bed height, the equations used to determine the concentration and velocity profiles in the sheet flow layer cannot be valid anymore. So the model is limited to a certain spatial volumetric concentration. If theoretically there is a very sliding small bed with a sheet flow layer on top, one can doubt whether such a situation can exist.

    The model has been tested with a number of equivalent bed roughness equations, but it would also be interesting to test the model with different erosion rate equations.

    Resuming one can say that the Matousek (2009), (2011) and Matousek & Krupicka (2010), (2011) model is based on a new very interesting concept. The user should however be aware of the empirical basis and some limitations of the model.

    The new shear stress based concentration profile model is very promising, compared to the traditional turbulent diffusion models for stratified flows. For 100% suspended load one can already discuss whether turbulent diffusion models are appropriate for particles that to not follows the motions of eddies, in stratified flows the particles in contact in a stationary or sliding bed have no relation with turbulence anymore. One has to distinguish here between fluid dynamics (suspended load) and soil mechanics (stratified loads). The new shear stress based model is doing exactly this. Once the mixture hydraulic gradient is known from an appropriate model/method, the height of a stationary or sliding bed and the height of the transport layer (shear layer or sheet flow layer) can be determined, both not following the turbulent diffusion model. Above the transport layer, the concentration profile can be determined with a turbulent diffusion model.

    One of the main issues is that the Richardson & Zaki (1954) hindered settling equation is based on the spatial volumetric concentration Cvs and not on the relative spatial volumetric concentration Cvr=Cvs/Cvb.

    \[\ \frac{\mathrm{v}_{\mathrm{t h}}}{\mathrm{v}_{\mathrm{t}}}=\left(1-\mathrm{C}_{\mathrm{v s}}\right)^{\beta}\]

    So even when the spatial volumetric concentration reaches a concentration where a bed with maximum porosity occurs, for sand at about Cvs=50%, still a hindered settling velocity is determined, while in reality this hindered settling velocity will be zero. Normal sands will have a porosity of about 40%, so Cvb=60%. A fixed bed may have a porosity of 40%, but a sliding bed will have a higher porosity in between 40% and 50%. The porosities mentioned here depend on the type of sand, but are mentioned to give a feeling of the order of magnitude. Since the Richardson & Zaki (1954) equation is based on small concentrations it is better to use a modified equation based on the relative concentration, for example:

    \[\ \frac{\mathrm{v}_{\mathrm{t h}}}{\mathrm{v}_{\mathrm{t}}}=\left(1-\mathrm{C}_{\mathrm{v r}}\right)^{\beta^{\prime}}\]

    Of course the power of this equation will be different from the original equation. An equation that may even work better is:

    \[\ \frac{\mathrm{v}_{\mathrm{t h}}}{\mathrm{v}_{\mathrm{t}}}=\mathrm{e}^{-\beta \cdot \mathrm{C}_{\mathrm{vr}}^{1.25}} \cdot\left(1-\mathrm{C}_{\mathrm{v r}}^{2}\right)\]

    For small concentrations this equation gives the same result as the original equation, but for concentrations approaching the bed concentration, this equation approaches a zero settling velocity. This would describe the bed behavior much better. So for small concentrations this equation describes hindered settling, while for large relative concentrations approaching 1, the behavior is more close to consolidation. The power β in this equation is equal to the original power β.

    6.24.6 Nomenclature Matousek Model


    Cross section pipe



    Cross section above bed



    Wall associated cross section upper layer



    Bed associated cross section upper layer



    Cross section bed



    Cross section sheet flow layer



    Volumetric spatial bed concentration



    Spatial volumetric concentration



    Delivered volumetric concentration in cross section 1, upper layer



    Delivered volumetric concentration in cross section 2, lower layer



    Delivered (transport) volumetric concentration



    Particle diameter



    Particle diameter with 50% passing



    Pipe diameter



    Relative excess hydraulic gradient



    Force on bed due to sliding friction



    Force on bed due to pressure and shear stress on top of the bed



    Gravitational constant 9.81 m/s2



    Thickness sheet flow layer



    Hydraulic gradient liquid



    Hydraulic gradient mixture



    Bed roughness



    Length of pipe section



    Limit Deposit Velocity



    Limit of Stationary Deposit Velocity



    Power of velocity distribution in sheet flow layer



    Circumference pipe



    Circumference pipe above bed


    O2 ​​​​​​​

    Circumference pipe in bed



    Width of bed



    Total solids flow rate


    Qs1 ​​​​​​​

    Solids flow rate in the upper layer


    Qs2 ​​​​​​​

    Solids flow rate in the lower layer, the bed



    Reynolds number



    Particle Reynolds number


    Relative submerged density



    Hydraulic radius upper layer



    Wall associated hydraulic radius upper layer



    Bed associated hydraulic radius upper layer



    Hydraulic radius lower layer



    Friction velocity



    Terminal settling velocity



    Line speed


    v1 ​​​​​​​

    Cross section averaged velocity above bed


    v2 ​​​​​​​

    Cross section averaged velocity bed



    Height of bed


    α ​​​​​​​

    MPM coefficient


    β ​​​​​​​

    MPM coefficient



    Bed angle



    Pipe wall roughness



    Internal friction angle


    \(\ \gamma\)

    Proportionality constant sheet flow layer



    Density carrier liquid



    Density solids


    \(\ \kappa\)

    Von Karman constant, about 0.4



    Shields parameter



    Critical Shields parameter



    Darcy-Weisbach friction factor liquid-pipe wall



    Darcy-Weisbach friction factor with pipe wall



    Darcy-Weisbach friction factor with pipe wall, liquid in bed



    Darcy-Weisbach friction factor on the bed


    \(\ v_{\mathrm{l}}\)

    Kinematic viscosity


    \(\ \tau_{\mathrm{l}}\)

    Shear stress liquid-pipe wall


    \(\ \mathrm{\tau}_{1,\mathrm{l}}\)

    Shear stress liquid-pipe wall above bed


    \(\ \tau_{12,\mathrm{l}}\)

    Shear stress bed-liquid


    \(\ \tau_{2,\mathrm{l}}\)

    Shear stress liquid-pipe in bed


    \(\ \tau_{2, \mathrm{sf}}\)

    Shear stress from sliding friction



    Sliding friction coefficient



    Slip ratio


    This page titled 6.24: The Matousek (2010), (2011) Model is shared under a CC BY-NC-SA 4.0 license and was authored, remixed, and/or curated by Sape A. Miedema (TU Delft Open Textbooks) via source content that was edited to the style and standards of the LibreTexts platform; a detailed edit history is available upon request.