∂q           ∂d            _____  +  _____  =  0   ∂x           ∂t

(1)

Continuidad de sedimentos:

∂gs         ∂                                       ∂z ____  +  ____ (Cs d)  +  (1 - p) ρs g ____  =  0  ∂x          ∂t                                       ∂t

(2)

Ecuación de movimiento:

1       ∂u           u       ∂u            ∂h             τb ____  _____  +   ____  _____  +   _____  +   _______  =  0   g       ∂t            g       ∂x            ∂x           ρw gd

(3)

en las cuales: u = velocidad media; h = elevación de la superficie del agua; z = elevación del fondo del canal; d = profundidad de flujo = h - z; q = descarga, o caudal, por unidad de ancho = ud; g_s = función de transporte de material de fondo, por unidad de ancho; p = porosidad del material depositado en el lecho; g = aceleración de la gravedad ρ_w = densidad del agua; ρ_s = densidad de sólidos; τ_b = esfuerzo cortante de fondo; x = variable de espacio, y t = variable de tiempo.

Se necesitan dos ecuaciones complementarias que describan la fricción del fondo y el transporte de material de fondo. La fricción de fondo se puede expresar mediante una relación de tipo Chézy, de la siguiente forma:

τb = ƒρw u 2

(4)

en la cual ƒ = factor de fricción adimensional, definido como sigue:

g ƒ =  _____          C 2

(5)

y C = coeficiente de Chézy. Para flujo uniforme, el esfuerzo cortante de fondo es:

τbo = ƒρw uo2

(6)

De igual manera:

τbo = ρw g do So

(7)

en la cual S_o = pendiente de fondo. Por lo tanto, a partir de las Ecs. 6 y 7 se obtiene:

So = ƒFo2

(8)

en la cual F_o = número de Froude del flujo uniforme, definido como sigue:

uo Fo =  __________            (g do) 1/2

(9)

La velocidad de transporte del material de fondo se modela mediante una relación de tipo Colby (1) de la siguiente forma:

gs = k um

(10)

en la cual k y m son parámetros empíricos de transporte de material de fondo.

3.  ECUACIONES SIMPLIFICADAS

En el flujo subcrítico de régimen bajo (F_o ≤ 0.6), las ondas transitorias del lecho se propagan aguas abajo con una celeridad que es varias órdenes de magnitud menor que la de las ondas de la superficie. Esta diferencia de escala de tiempo permite la solución para las ondas del lecho asumiendo una descarga constante dentro de un intervalo de tiempo, reduciendo así el problema a dos ecuaciones (continuidad de sedimentos y ecuación de movimiento) en dos incógnitas (nivel de la superficie del agua y nivel del lecho de la corriente). Además, para F_o ≤ 0.6, el término de concentración espacial en la Ec. 2 es pequeño en comparación con los términos restantes (4); por lo tanto, puede despreciarse en una primer aproximación. Para facilitar el análisis, y en parte sobre la base de un razonamiento de orden de magnitud, el término de aceleración local en la Ec. 3 también puede despreciarse.

Sustituyendo la Ec. 10 en la Ec. 2, y expresando u en términos de q y d:

∂                                           ∂z k qm _____ ( d -m ) +  ( 1 - p ) ρs g _____  =  0           ∂x                                           ∂t

(11)

Operando en la Ec. 11:

m k qm        ∂d                                ∂z ( - _________ ) ______  +  ( 1 - p ) ρs g _____  =  0        d m+1          ∂x                                 ∂t

(12)

Sustituyendo la Ec. 10 en la Ec. 12:

m gs         ∂d                                ∂z ( - _______ ) ______  +  ( 1 - p ) ρs g _____  =  0          d            ∂x                                 ∂t

(13)

Dividiendo por u ( 1 - p ) ρ_s g:

m cb                ∂d               1          ∂z [ - ______________ ] ______  +  ( _____ ) ______  =  0                     ρs            ∂x               u          ∂t       ( 1 - p ) ____                     ρw

(14)

en la cual c_b = concentración de material del lecho en partes por unidad, definida como sigue:

gs cb  =  ___________             ρw g u d

(15)

En la Ecuación 14, expresando u en términos de q y d, y d en términos de h y z:

m cb                ∂h                    m cb                ∂z             h - z          ∂z [ - ______________ ] ______  +  [ ______________ ] ______  + ( _______ ) ______ =  0                     ρs            ∂x                             ρs          ∂x                q            ∂t       ( 1 - p ) ____                             ( 1 - p ) ____                     ρw                                            ρw

(16)

Sustituyendo la Ecuación 4 en la Ecuación 3:

u      ∂u        ∂h         ƒu2 ____ _____ + _____ + _____ =  0   g      ∂x        ∂x         gd

(17)

En la Ecuación 17, expresando u en términos de q y d, y d en términos de h y z:

q 2        ∂h             q 2        ∂z         ƒq 2 ( 1 - ______ ) _____ + ( ______ ) _____ + ______ =  0          g d 3       ∂x           g d 3       ∂x        g d 3

(18)

4.  ANÁLISIS DE ESTABILIDAD LINEAL

Las Ecuaciones 16 y 18 deben satisfacer tanto el flujo de equilibrio para el cual h = h_o, z = z_o, c_b = c_bo, y d = d_o, y el flujo transitorio, para el cual h = h_o + h^', z = z_o + z^', c_b = c_bo + c^'_b, y d = d_o + d^'. El superíndice ^' representa una pequeña perturbación del flujo de equilibrio. Por lo tanto, todos los términos cuadráticos en los componentes perturbados pueden despreciarse debido a un razonamiento de orden de magnitud.

La sustitución de las variables perturbadas en la Ec. 16 conduce a la siguiente ecuación:

∂h'              ∂z'          ho - zo        ∂z' -Ao _____ + Ao ______ + __________ ______ =  0         ∂x               ∂x              q             ∂t

(19)

en la cual:

m cbo Ao  = _______________                          ρs             (1 - p) ______                          ρw

(20)

De las Ecuaciones 9, 10 y 15:

k Fo 2 uo m - 3 cbo  =  ________________                       ρw

(21)

Por lo tanto:

Ao = Φ Fo 2

(22)

en la cual Φ = parámetro compuesto de transporte de material de lecho, definido como sigue:

m k uo m - 3 Φ = ______________            (1 - p) ρs

(23)

Sustituyendo la Ec. 22 en la Ec. 19:

∂h'                    ∂z'         ho - zo      ∂z' -Φ Fo2 _____  + Φ Fo2 _____  + _________ _____  = 0              ∂x                     ∂x             q           ∂t

(24)

La sustitución de las variables perturbadas en la Ec. 18 conduce a lo siguiente:

∂h'                ∂z'                      h' - z'       (1 - Fo 2) _____  + Fo2 _____ - 3ƒFo2 ( _________ ) = 0                  ∂x                ∂x                      ho - zo

(25)

en la cual:

q 2 Fo 2  =  _______               g do3

(26)

5.  SOLUCION ANALÍTICA

La solución del sistema de Ecs. 24 y 25 se postula en forma exponencial como sigue:

h' ____ = h* exp [i (σ* x* - β* t*)]   do

(27)

y

z' ____ = z* exp [i (σ* x* - β t*)]   do

(28)

en las cuales:

2π σ* = ( _____ ) Lo              L

(29)

2π        Lo β*R =  ( _____ ) _____                 T         uo

(30)

β*I = factor de atenuación (o amplificación)

(31)

x x* = ____          Lo

(32)

uo t* = t _____            Lo

(33)

y

do Lo = _____           So

(34)

en las cuales L = longitud de onda; T = período de onda; y h_* y z_* son funciones adimensionales de amplitud.

La celeridad adimensional, c_*, de una perturbación sinusoidal se define como sigue (6):

L         ___          T         β*R c* = ____ = ______          uo         σ*

(35)

El decremento logarítmico, δ, se define de la siguiente manera (6):

β*I δ = 2π ________              | β*R |

(36)

La sustitución de las Ecs. 27 y 28 en las Ecs. 24 y 25 dan como resultado un sistema homogéneo de ecuaciones lineales en las incógnitas h_* y z_*. La condición no trivial para el determinante de la matriz de coeficientes conduce a una celeridad adimensional normalizada de transporte, c_*aΦ, expresada como sigue:

c*a            Fo2(1 - Fo2)σ*2 c*aΦ = _____ = ____________________               Φ            (1 - Fo2)σ*2 + 9

(37)

y un decremento logarítmico, δ_a, expresado como sigue:

6π δa = - ____________             |1 - Fo2|σ*

(38)

Las Figuras 1 y 2 muestran la celeridad adimensional normalizada por transporte, c_*aΦ, y el decremento logarítmico, δ_a, en función del número de Froude del flujo de equilibrio, F_o, y el número de onda adimensional, σ_*, en flujo subcrítico.

Fig. 1  Celeridad analítica adimensional normalizada por transporte c*aΦ versus número de onda adimensional σ*, en flujo subcrítico.

Fig. 2  Decremento logarítmico de la solución analítica δa versus número de onda adimensional σ*, en flujo subcrítico.

---

6.  SOLUCIÓN NUMÉRICA

El esquema implícito de cuatro puntos de Preissmann (3) se utiliza para discretizar las ecuaciones de gobierno. Las ecuaciones en diferencias finitas tienen la siguiente forma (Fig.3):

Fig. 3  Bosquejo de definición para el esquema implícito de cuatro puntos.

θ                                  1 - θ ƒ(x,t) = ___ (ƒj + 1n + 1 + ƒjn + 1) + ______ (ƒj + 1n + ƒjn)               2                                     2

(39)

∂                   1 ____ ƒ(x,t) ≈ _____ (ƒj + 1n + 1 - ƒj + 1n + ƒjn + 1 - ƒjn)  ∂t                  2∆t

(40)

∂                   θ                               1 - θ ____ ƒ(x,t) ≈ _____ (ƒj + 1n - ƒjn + 1) + ______ (ƒj + 1n - ƒjn)  ∂x                 ∆x                                ∆x

(41)

en las cuales ƒ(x,t) = variable dependiente; Δx = intervalo de espacio; Δt = intervalo de tiempo; y θ = factor de ponderación.

Llamando η = h^' y ζ = z^', las Ecs. 24 y 25 se convierten a:

∂η                  ∂ζ        do      ∂ζ -Φ Fo2 _____ + Φ Fo2 ____ + _____ _____ = 0              ∂x                  ∂x         q       ∂t

(42)

∂η               ∂ζ                   η - ζ     (1 - Fo2) _____ + Fo2 ____ - 3ƒFo2 (_______) = 0                 ∂x               ∂x                     do

(43)

Sustituyendo las Ecs. 39-41 en 42 y 43 se obtiene:

θ                               1 - θ                                      θ -ΦFo2[_____ (ηj+1n+1 - ηjn+1) + ______ (ηj+1n - ηjn)] + ΦFo2[_____ (ζj+1n+1 - ζjn+1)              Δx                                Δx                                      Δx      1 - θ                            do + _______ (ζj+1n - ζjn)] + _______ [ζj+1n+1 - ζj+1n + ζjn+1 - ζjn] = 0        Δx                           2qΔt

(44)

θ                               1 - θ                                   θ (1 - Fo2)[_____ (ηj+1n+1 - ηjn+1) + ______ (ηj+1n - ηjn)] + Fo2[_____ (ζj+1n+1 - ζjn+1)                  Δx                                Δx                                   Δx      1 - θ                       -3ƒFo2       θ + _______ (ζj+1n - ζjn)] _______ {[ ____ (ηj+1n+1 + ηjn+1)        Δx                          do            2      1 - θ                               θ                                1 - θ + _______ (ηj+1n + ηjn)] - [ _____ (ζj+1n+1 + ζjn+1) + _______ (ζj+1n + ζjn) ]} = 0         2                                  2                                   2

(45)

La sustitución de las Ecs. 27 y 28 en las Ecs. 44 y 45 conduce a lo siguiente:

Δt*                                                   Δt*                     E h*[ -ΦFo2 (θE + 1)( _____ ) i tanα]+ z*[ ΦFo2 (θE + 1)( _____ ) i tanα + ____ ] = 0                                   Δx*                                                  Δx*                     2

(46)

Δt*                   3 h*[ (1 - Fo2) (θE + 1)( _____ )i tanα - ____(θE + 1)Δt*]                                       Δx*                  2                                  Δt*                     3 + z*[ Fo2 (θE + 1)( _____ ) i tanα + ____(θE + 1)Δt*] = 0                                  Δx*                    2

(47)

en las cuales:

E = e -iβ*Δt* - 1

(48)

Δx Δx* = _____             Lo

(49)

uo Δt* = Δt _____                 Lo

(50)

y

σ*Δx* α = ________             2

(51)

Las Ecuaciones 46 y 47 conducen a un sistema homogéneo de ecuaciones lineales con las incógnitas h_* y z_*. Para que la solución no sea trivial, el determinante de la matriz de coeficientes debe anularse. Esto conduce a:

1        c*a                                        3        c*a E [ Φθσ*Fo2 tan2α - ____ ( _____ ) σ* (1 - Fo2) i tanα + ____ ( _____ ) α ]                                    2         C                                         2         C + Φσ*Fo2 tan2α = 0

(52)

en la cual C = número de Courant del lecho, definido como sigue:

L          ___           T               Δt C = ______ = ca _____          Δx              Δx         ____          Δt

(53)

Dado que:

c*a c*aΦ = _____               Φ

(54)

la Ecuación 52 se reduce a:

1       c*aΦ                                       3      c*aΦ E [ Φσ*Fo2 tan2α - ___ ( ______ ) σ* (1 - Fo2) i tanα + ___ ( ______ ) α ] + σ*Fo2 tan2α = 0                                 2         C                                         2        C

(55)

Resolviendo la Ec. 55 para E:

M + i N exp (-iβ*Δt*)  =  1 - ________                                     P

(56)

en la cual:

3       c*aΦ M = θ (σ*Fo2 tan2α)2 + ___ ( ______ ) (σ*Fo2 α tan2α)                                        2        C

(57)

1       c*aΦ N = ___ ( ______ ) σ*2 (1 - Fo2)Fo2 tan3α         2         C

(58)

3       c*aΦ                  1       c*aΦ P = [ θσ*Fo2 tan2α + ___ ( ______ )α ]2 + [ ___ ( ______ )σ* (1 - Fo2) tanα]2                                     2         C                     2         C

(59)

En la Ecuación 56, expresando β_* en sus componentes reales e imaginarios:

β* = β*R + i β*I

(60)

Se obtiene:

N tan (β*R Δt*) = _______                          P - M

(61)

y

[(P - M) 2 + N 2]1/2 exp(β*I Δt*) = ___________________                                     P

(62)

La celeridad adimensional de la solución numérica, c_*n, se define como sigue:

N                              tan-1( ________ )             β*R                        P - M c*n = ________ = ___________________               σ*                    σ* Δt*

(63)

Dado que:

2α σ* Δt* = ______                 c*a                _____                  C

(64)

y definiendo:

c*n c*nΦ = _____               Φ

(65)

la Ecuación 63 se reduce a:

c*aΦ                     N             ( ______ ) tan-1 ( _______ )                   C                    P - M c*nΦ = ___________________________                             2α

(66)

El decremento logarítmico de la solución numérica, δ_n, es definido como sigue:

[(P - M)2 + N 2]1/2                 β*I               ln { ____________________ }                                                         P δn = 2π _______ = 2π _____________________________               |β*R |                                    N                                          |tan-1 ( _______ )|                                                        P - M

(67)

Las Ecuaciones 66 y 67 permiten el cálculo de la celeridad de propagación, c_*nΦ, y decremento logarítmico, δ_n, del esquema numérico implícito de cuatro puntos de las Ecs. 2 y 3, en función de: (1) Número de Froude del flujo de equilibrio, F_o = u_o / (gd_o)^1/2; (2) número de onda adimensional, σ_*, del componente transitorio del movimiento, σ_* = (2π /L)L_o; (3) parámetro de resolución espacial, α = π / (L / Δx); (4) número de Courant del lecho móvil, C = c_a / (Δx / Δt); y (5) factor de ponderación, θ.

---

7.  ANÁLISIS DE CONVERGENCIA

El análisis de convergencia se basa en las siguientes ralaciones:

R1 = exp (δn - δa)

(68)

c*nΦ R2 = _______           c*aΦ

(69)

El valor R₁ es la relación de atenuación, y R₂ es la relación de translación (fase). Para (δ_n - δ_a) > 0, R₁ > 1, lo cual conduce a amplificación numérica; para (δ_n - δ_a) < 0, R₁ < 1, lo cual conduce a una atenuación numérica. Para R₂ > 1, la translación numérica es mayor que la translación física; para R₂ < 1, la translación numérica es menor que la translación física. Para R₁ = R₂ = 1, existe una coincidencia exacta entre la soluciones analíticas y numéricas.

Las Figuras 4 y 5 muestran la relación de atenuación, R₁, como una función de F_o, σ_*, α, C y θ. En la Fig. 4, α se mantiene constante e igual a π / 40; en la Fig. 5, σ_* se mantiene constante e igual a 10. Las Fig. 4 y 5 muestran que al aumentar σ_* y disminuir F_o, α, y C da como resultado una mejor convergencia de la amplitud de onda, es decir, R₁ → 1. Un valor de θ = 0.5 es una condición necesaria para la estabilidad (R₁ < 1), pero no es suficiente para el caso general. La inestabilidad (R₁ > 1) puede observarse para ciertas combinaciones de F_o, σ_*, y α, aún cuando θ > 0.5. Interesantemente, este comportamiento se ha observado en el funcionamiento del modelo real (2, 4, 7).

Fig. 4  Relación de atenuación R1 en función de Fo, σ*, C, y θ; α = π/40.

Fig. 5  Relación de atenuación R1 en función de Fo, α, C, y θ; σ* = 10.

Las Figuras 6 y 7 muestran la relación de translación, R₂, en función de F_o, σ_*, α, C y θ. En la Fig. 6, α se mantiene constante e igual a π / 40; en la Fig. 7, σ_* se mantiene constante e igual a 10. Nuevamente, las Figs. 6 y 7 muestran que al aumentar σ_* y disminuir F_o, α, y C, resulta en una mejor convergencia de la translación de onda, es decir, R₂ → 1.

Fig. 6  Relación de translación R2 en función de Fo, σ*, C, y θ; α = π/40.

Fig. 7  Relación de translación R2 en función de Fo, α, C, y θ; σ* = 10.

Los resultados de las Figs. 4-7 permiten obtener las siguientes conclusiones: (1) La convergencia mejora para un valor suficientemente alto de σ_*, y para valores bajos de F_o, α, y C; y (2) para un valor de θ en el rango medio, entre 0.5 y 1, satisfacerá en general los requisitos de estabilidad y convergencia. Con base en el análisis realizado, se recomienda el siguiente rango de parámetros: σ_* ≥ 10; F_o ≤ 0.6; α ≤ π / 40; C ≤ 2; y 0.6 ≤ θ ≤ 0.8.

En la práctica, σ y F_o están determinadas por las características del problema siendo considerado. Los valores recomendados para α, C, y θ deben tomarse solo como una sugerencia. La experimentación numérica es necesaria para evaluar la sensibilidad de los resultados del modelo a los parámetros de convergencia. Los experimentos numéricos tienen como objetivo lograr una convergencia satisfactoria (es decir, minimizar la dispersión numérica) mientras se mantiene la estabilidad.

---

8.  RESUMEN Y CONCLUSIONES

Se presenta un tratamiento teórico detallado de la convergencia del modelo numérico implícito transitorio de lecho móvil. Se calcula la celeridad adimensional y el decremento logarítmico de las soluciones analíticas y numéricas. La convergencia se prueba estableciendo la relación de celeridades (convergencia de translación) y de las relaciones de atenuación (convergencia de atenuación). Se identifican dos parámetros físicos y tres numéricos, a saber: el número de Froude del flujo de equilibrio F_o, el número de onda adimensional del flujo transitorio σ_* , la resolución espacial α, el número de Courant C, y el factor de ponderación θ.

---

AGRADECIMIENTOS

La investigación presentada en este documento fue apoyada en parte por la Estación Experimental de la Universidad Estatal de Colorado. La visita de Horst Indlekofer a la Universidad Estatal de Colorado fue posible gracias a una beca de la OTAN.

---

APÉNDICE Ι. BIBLIOGRAFÍA

1. Colby, B, R. 1964. "Discharge of Sands and Mean Velocity Relationships in Sand-Bed Streams," Professinal Paper 462-A, United States Geological Survey, Washington, D.C.

2. Cunge, J. A., y N. Perdreau. 1973. "Mobile Bed Fluvial Mathematical Models," La Houille Balanche, Grenoble, France, No. 7, 561-580.

3. Liggett, J. A., y J. A. Cunge. 1975. "Numerica Methods of Solution of the Unsteady Flow Equations," Unsteady Flow in Open Channels, K. Mahmood and V. Yevjevich, eds., Water Resources Publications, Fort Collins, Colo., 89-182.

4. Mahmood, K., and V. M. Ponce. 1976. "Mathematical Modeling of Sedimentation Transients in Sand-Bed Channels," Report CER75-76-KM-VMP28, Engineering Research Center, Colorado State University, Fort Collins, Colo., Apr.

5. O'Brien, G. G., M. A. Hyman, and S. Kaplan. 1951. "A study of the Numerical Solution of Partial Differential Equations," Journal of Mathematics and Physics, Vol. 29, No. 4, 223-251.

6. Ponce, V. M., y D. B. Simons. 1977. "Shallow Wave Propagation in Open Channel Flow," Journal of the Hydraulics Division, ASCE, Vol. 103, No. HY12. Proc. Paper 13392, Dec., 1461-1476.

7. Ponce, V. M., J. Lopez Garcia, y D. B. Simons. 1979. "Modeling Alluvial Channel Bed Transients," Journal of the Hydraulics Division, ASCE, Vol. 105, No. HY3. Proc. Paper 14456, Mar., 245-256.

8. Roache, P. J. 1972. Computational Fluid Dynamics, Hermosa Publishers, Albuquerque, N. M.

---

APÉNDICE ΙΙ. SIMBOLOGÍA

En este artículo se utilizan los siguientes símbolos:

A_o = parámetro, Ec. 20;

C = coeficiente de Chézy, y número Courant de lecho móvil;

c_a = celeridad analítica;

c_b = concentración de material del lecho, Ec. 15;

c_n = celeridad numérica;

c_* = celeridad adimensional, Ec. 35;

c_*a = celeridad analítica adimensional;

c_*n = celeridad numérica adimensional;

c_*aΦ = celeridad analítica adimensional normalizada de transporte, Ec. 37;

c_*nΦ = velocidad numérica adimensional normalizada de transporte, Ec. 66;

d = profundidad de flujo;

E = parámetro, Ec. 48;

F = número de Froude;

ƒ = factor de fricción adimensional, Ec. 5 y variable dependiente;

g = aceleración de la gravedad;

g_s = tasa de transporte de material del lecho de ancho unitario;

h = elevación de la superficie del agua;

j = índice de discretización espacial;

k = parámetro empírico de transporte de material del lecho, Ec. 10;

L = longitud de onda;

L_o = longitud del canal, Ec. 34;

M = parámetro, Ec. 57;

m = parámetro empírico de transporte de material del lecho, Ec. 10;

N = parámetro, Ec. 58;

n = índice de discretización de tiempo;

P = parámetro, Ec. 59;

p = porosidad del material del lecho;

q = descarga de agua, por unidad de ancho;

R₁ = relación de atenuación, Ec. 68;

R₂ = relación de translación (fase), Ec. 69;

S_o = pendiente de fondo;

T = período de la onda;

t = variable temporal;

u = velocidad media;

x = variable espacial;

z = elevación de fondo;

α = parámetro de resolución espacial, Ec. 51;

β_*I = factor de atenuación (o amplificación);

β_*R = frecuencia adimensional, Ec. 30;

Δt = intervalo de tiempo;

Δx = intervalo de espacio;

δ = decremento logarítmico, Ec. 36;

δ_a = decremento logarítmico de la solución analítica, Ec. 38;

δ_n = decremento logarítmico de solución numérica, Ec. 67;

ζ = elevación de la perturbación de fondo;

η = perturbación de la elevación de la superficie del agua;

θ = factor de ponderación del esquema implícito;

ρ_s = densidad de sólidos;

ρ_w = densidad del agua;

σ_* = número de onda adimensional, Ec. 29;

τ_b = esfuerzo cortante de fondo;

Φ = parámetro compuesto de transporte de material de lecho, Ec. 23;

Subíndices

o = flujo de equilibrio;

* = variable adimensional; y

Superíndice

' = variable perturbada.

---