Full text
An´ alisis de modelos de fricci´ on en flujos de superficie libre Pedro Mart´ın Navarro M´aster en Mec´anica Aplicada Programa Oficial de Posgrado en Ingenier´ıa Mec´anica y de Materiales Noviembre 2010 Director: Javier Murillo Castarlenas Centro Polit´ecnico Superior. Universidad de Zaragoza
. 2
An´ alisis de modelos de fricci´ on en flujos de superficie libre Resumen La modelizaci´on de las p´erdidas por fricci´on que se producen en un flujo turbulento bidimensional de superficie libre se ha basado tradicionalmente en la extensi´on simplista de ideas desarrolladas para flujos unidimensionales en condiciones de equilibrio din´amico. Es preciso establecer las condiciones de aplicabilidad de estos modelos e investigar la posible superioridad de otras formulaciones. Para poder generar modelos m´as ambiciosos que los t´ıpicamente utilizados en la hidr´aulica fluvial, es ´util plantear adecuadamente el problema de la descripci´on del esfuerzo en el fondo. Por este motivo, es necesario recopilar la informaci´on existente acerca de modelos de p´erdidas por fricci´on en el contexto de flujos de superficie libre promediados en la vertical. La teor´ıa de capa l´ımite, que nos permite describir la variaci´on de la velocidad en la direcci´on vertical, es el instrumento adecuado. Su interpretaci´on permite comprender otras aproximaciones, como las basadas en perfiles potenciales de velocidad. El an´alisis La comprensi´on de los diferentes modelos de fricci´on puede facilitarse a trav´es de la simulaci´on num´erica que, adem´as, es el campo de aplicaci´on que se persigue. Comparar la capacidad de predicci´on del comportamiento de flujos bidimensionales estacionarios caracterizados experimentalmente es un paso fundamental, en torno al cual gira este trabajo. 3
. 4
´ Indice 1. Introducci´on 7 2. Ecuaciones gobernantes del flujo de superficie libre 8 3. Capa l´ımite y esfuerzo en el fondo 12 3.1. Aplicaci´on a fondo rugoso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2. Esfuerzo en lechos fluviales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3. La aproximaci´on de Manning-Strickler . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.4. La aproximaci´on de Burguete . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4. Modelos de fricci´on en canales compuestos 21 5. An´alisis de los datos experimentales 23 6. Resultados mediante simulaci´on num´erica 2D 34 7. Conclusiones y trabajo futuro 47 8. Bibliograf´ıa 48 9. Anexo: El esquema num´erico 52 5
. 6
1. Introducci´on Han sido muchos los que han aportado con sus investigaciones al modelado de la fricci´on en el lecho de un canal o r´ıo, con distintos enfoques [2, 7, 12, 16, 22]. Los primeros que investigaron y relacionaron par´ametros como la pendiente de carga, y la velocidad en la secci´on propusieron modelos unidimensionales con un solo par´ametro dependiente del material y la geometr´ıa del conducto y hoy en d´ıa este punto de vista sigue siendo ampliamente usado en el campo de la hidr´aulica. Parece razonable que estas primeras investigaciones, tanto en tuber´ıas como en conductos abiertos, empezaran modelando el sistema hidr´aulico de forma unidimensional, es decir: las variables f´ısicas son forzadas a tener una sola componente en la direcci´on del flujo, lo que es l´ogico si se tiene en cuenta que se est´a hablando de canales y por tanto de flujos encauzados y de secci´on conocida, con una geometr´ıa que se repite a lo largo de un eje. Nos expresamos por tanto en t´erminos de fricci´on por ´area de secci´on y velocidad promediada en la secci´on. Actualmente, una rama importante de estudio es la predicci´on, a trav´es de la simulaci´on num´erica, de avenidas de inundaci´on producidas por r´ıos en cuencas hidrol´ogicas. Esto es equivalente a tener canales con al menos dos fondos de alturas distintas, ya que se quiere modelar el cauce de un r´ıo con llanuras de inundaci´on a ambos lados. El flujo que circula por el cauce acelerar´a el de las llanuras, y ´este decelerar´a aqu´el, lo que nos advierte de una transferencia de cantidad de movimiento en el eje perpendicular al de avance del flujo [20]. Y por tanto, de la necesidad de un estudio bidimensional si se quiere recoger lo que est´a ocurriendo a lo largo de las diferentes secciones. Para poder generar modelos m´as ambiciosos que los t´ıpicamente utilizados en la hidr´aulica fluvial es necesario plantear adecuadamente el problema de la descripci´on del esfuerzo en el fondo. Por este motivo, en la primera parte del trabajo se estudian las propiedades de la capa l´ımite que nos permiten describir la variaci´on de la velocidad en la direcci´on vertical, acotada por la altura de la l´amina de agua. La interpretaci´on de los resultados de la capa l´ımite junto con nuevas aproximaciones basadas en perfiles potenciales de velocidad, sirve en la segunda parte del trabajo para evaluar la idoneidad de ampliar los modelos cl´asicos con la m´ınima informaci´on posible. 7
2. Ecuaciones gobernantes del flujo de superficie libre En el modelo de flujo llamado de aguas poco profundas el hecho esencial es que el grosor h de la capa de fluido es peque˜no comparado con la escala longitudinal horizontal t´ıpica L(Fig. 1). Este es el punto en com´un que guardan todas las aplicaciones basadas en el modelo. Las ecuaciones que describen el flujo tridimensional de superficie libre [35] pueden ser obtenidas a partir de las ecuaciones de Navier-Stokes que expresan los principios f´ısicos de conservaci´on de la masa y conservaci´on del movimiento en las tres direcciones del espacio (1), (2). Si las variaciones de densidad son despreciables o nulas, no existen v´ınculos entre la ecuaci´on de conservaci´on de la energ´ıa y las de la masa y movimiento y en este caso son suficientes las dos ´ultimas ecuaciones para conocer la evoluci´on del flujo. Para entender el alcance de las hip´otesis asociadas hay que deducir las ecuaciones de aguas poco profundas a partir de la ecuaci´on de conservaci´on de la masa ∂u ∂x +∂v ∂y +∂w ∂z = 0 (1) donde u, v, w son las tres componentes de la velocidad del flujo en las direcciones x, y, z respectivamente, y de las ecuaciones de Navier-Stokes para flujo incompresible, que en forma diferencial conservativa se escriben como: ∂ρu ∂t +∂ρu2 ∂x +∂ρuv ∂y +∂ρuw ∂z =−g∂p ∂x +∂τxx ∂x +∂τxy ∂y +∂τxz ∂z ∂ρv ∂t +∂ρuv ∂x +∂ρv2 ∂y +∂ρvw ∂z =−g∂p ∂y +∂τxy ∂x +∂τyy ∂y +∂τyz ∂z ∂ρw ∂t +∂ρuw ∂x +∂ρvw ∂y +∂ρw2 ∂z =−ρg −∂p ∂z +∂τxz ∂x +∂τyz ∂y +∂τzz ∂z (2) con ρla densidad del fluido, gel valor de la aceleraci´on en la vertical, τij las componentes del tensor de esfuerzos viscosos y pla presi´on. Todos los pasos a partir de estas ecuaciones deben formularse cuidadosamente, de tal forma que quede claro d´onde las aproximaciones son necesarias. L h xy z Figura 1: Dimensiones t´ıpicas del problema. Para fijar las soluciones de las ecuaciones diferenciales (1) y (2) es preciso definir las condiciones de contorno en las fronteras. En este caso existen dos zonas: la entrefase fluido-s´olido (fondo) que es fija y la entrefase fluido-fluido (superficie libre) que puede cambiar continuamente. 8
En el fondo la velocidad normal a la superficie s´olida debe ser cero (fondo s´olido, impermeable y fijo). vˆ nb=u∂zb ∂x +v∂zb ∂y −w= 0 (3) con ˆ nb= (∂zb/∂x, ∂zb/∂y, −1) el vector normal a la superficie l´ıquida en contacto con la s´olida hacia afuera en z=zb(x, y), donde zbes la cota del fondo medida desde un nivel de referencia horizontal (Fig. 2). Si el flujo es viscoso y el fondo fijo, la fuerza que act´ua es la de la viscosidad y por lo tanto las part´ıculas que se encuentran en contacto con el fondo est´an pegadas a ´el, por lo cual puede imponerse la condici´on de no deslizamiento que conduce a: u=v= 0 (4) en z=zb(x, y). Zb h H Figura 2: Perfil del cauce. En la superficie libre la velocidad normal relativa debe ser nula y representa que el fluido no se puede atravesar la superficie libre, situada en z=H(x, y, t), donde Hes el nivel de la superficie libre medido desde el nivel de referencia (Fig. 2). En la superficie libre se supone continuidad de esfuerzos. En este trabajo no van a ser considerados. Para especificar lo que se entiende por modelo simplificado de aguas poco profundas se tienen que considerar escalas t´ıpicas que permitan realizar un an´alisis de ´ordenes de magnitud. Las escalas longitudinales en la vertical de inter´es en este trabajo son: la profundidad de agua h el grosor δde la capa l´ımite en el fondo o en la superficie libre. Se generan por viscosidad y posiblemente se vean influenciadas por la rotaci´on de la tierra. Se considera una escala interna ya que es un par´ametro que viene determinado por el flujo. En r´ıos, lagos y costas, las capas l´ımite suelen ser m´as gruesas que la profundidad de agua, y entonces la capa de agua suele encontrarse incluida dentro de la capa l´ımite. En este caso, la profundidad hes la escala vertical relevante. 9
La ecuaci´on (34) se escribe generalmente utilizando la longitud z+y la velocidad u+definidas como z+=z u∗ νu+=u u∗ (35) con ν=µ/ρ y sustituyendo la constante multiplicadora por la constante del ingeniero h´ungaro von K´arm´an [36], que fue el precursor de esta t´ecnica. Introduciendo estos cambios, (34) se convierte en u+(z+) = 1 kln z++B(36) conocida como la ley logar´ıtmica v´alida en la regi´on intermedia y en la subcapa logar´ıtmica. La constante Bdepende del n´umero de Reynolds, y en el caso de flujos altamente turbulentos es una constante aproximadamente igual a 5.28. Aunque no puede derivarse f´acilmente de consideraciones dimensionales [19] el siguiente resultado para la subcapa viscosa puede derivarse: u+(z+) = (z+) (37) Las ecuaciones (36,37) son resultados plausibles para el caso de fondo liso o suave, aunque no totalmente consensuados en su rango de aplicaci´on [19]. 3.1. Aplicaci´on a fondo rugoso El flujo en r´ıos puede ser caracterizado por el esfuerzo cortante en el fondo, pero el fondo o contorno es raramente liso. Al contrario del caso suave, donde la ´unica forma de disipaci´on es la disipaci´on viscosa directa, en el caso rugoso el lecho experimenta fuerzas de arrastre en los elementos rugosos. Esta forma de p´erdida de energ´ıa cin´etica debe ser mantenida por las capas superiores del flujo. Esta es la caracter´ıstica m´as distintiva del flujo turbulento sobre lecho rugoso. Establecidas las diferencias, es necesario discernir que propiedades del caso suave pueden ser mantenidas. El grado de rugosidad del fondo puede caracterizarse a trav´es de la magnitud K+ s: K+ s=Ksu∗ ν(38) donde Kses el di´ametro de las part´ıculas fijas en una superficie, como en el caso del trabajo de Nikuradse, o una longitud din´amica equivalente. A este respecto, uno de los trabajos m´as relevantes fue el presentado por Raupach [28], donde se propone un mecanismo de producci´on de energ´ıa cin´etica turbulenta en flujo sobre paredes rugosas. Para ello, utiliz´o el concepto de subcapa rugosa, ´util para conceptualizar el balance de energ´ıa. La subcapa rugosa es la zona del flujo influenciada por los elementos rugosos salientes. La interacci´on entre las estelas y zonas de separaci´on detr´as de los elementos salientes rugosos es la caracter´ıstica fundamental de esta regi´on. Siguiendo trabajos previos [33], Raupach asumi´o que el movimiento turbulento es estrictamente independiente de la rugosidad (excepto en su influencia en la velocidad de fricci´on u∗), y apoy´andose en evidencias experimentales, demostr´o que la mayor diferencia entre un contorno liso y uno rugoso se confina dentro de la subcapa rugosa. Fuera de esta regi´on, las capas l´ımites lisas o rugosas tienen esencialmente la misma estructura. Las ecuaciones del perfil de velocidad en (36,37) se transforman en [8] 16
|ut| u∗ =1 kln z u∗E ν(39) donde utes la componente de velocidad paralela a la pared, zes la distancia normal a la pared y Ees un par´ametro de rugosidad. El par´ametro Edepende de la rugosidad adimensional K+ s E= 9,0 if K+ s≤5 1 0,11+0,033K+ sif 5 < K+ s<70 30 K+ sif K+ s≥70 (40) La evaluaci´on del par´ametro de rugosidad puede hacerse utilizando otras formulaciones [38]. Para paredes lisas Ks= 0. Para lechos rugosos de arena Ksse toma como el di´ametro medio del material del lecho, Ks=d50. Otras expresiones para Kspueden encontrase en [29]. 3.2. Esfuerzo en lechos fluviales La fricci´on en el fondo tiene un doble efecto en las ecuaciones promediadas en la vertical. Primero, produce una fuerza de fricci´on τbque se opone a la velocidad media del flujo. Y, segundo producen turbulencia. Ambos efectos se caracterizan por la velocidad de fricci´on u∗. Como en las ecuaciones promediadas no se resuelve en la vertical, es necesario relacionar la velocidad promedio ¯ ucon la velocidad cortante u∗, a trav´es de: |τb|=ρu2 ∗=ρcf|¯ u|2(41) Existen numerosas expresiones que permiten aproximar el coeficiente de fricci´on cf. La mayor´ıa de ellas suponen flujo completamente desarrollado y toda la columna de agua se aproxima como una ley logar´ıtmica de contorno. La rugosidad en r´ıos es apreciablemente alta, y puede asumirse que K+ s≥70, siendo su efecto importante, incluso a distancias comparativamente grandes del contorno fijo. En este caso, a partir de (39) se obtiene que |ut| u∗ =1 kln 30z Ks(42) Ahora la velocidad horizontal promedio puede obtenerse integrando (42) en la vertical |¯ u|=1 hZzs zb|ut|dz =u∗ kln 30h Ks−1 + Ks 30h(43) y considerando (41) se puede escribir cf=u2 ∗ |¯ u|2(44) obteniendo el siguiente valor de cf cf=1 kln 30h Ks−1 + Ks 30h−2 (45) 17
3.3. La aproximaci´on de Manning-Strickler Si se supone que la altura de la l´amina de agua es mucho mayor que la rugosidad del fondo, h >> Ksse obtiene la ecuaci´on de Keulegans: cf=1 kln 11h Ks−2 (46) que puede ser aproximada por c−1/2 f= 8,1h Ks1/6 (47) m´as conocida por la aproximaci´on de Manning-Strickler. Una forma alternativa de evaluar la fricci´on en el fondo, utilizada ampliamente en el hidra´ulica es la f´ormula de Manning [22] que utiliza el coeficiente de Manning, n, en vez de la altura de rugosidad, Ks. El coeficiente de fricci´on cfse escribe como: cf=gn2 h1/3(48) Comparando los coeficientes de fricci´on dados por la f´ormula de Manning (48) y la aproximaci´on de Manning-Strickler (47), aparece la siguiente relaci´on n≈K1/6 s 25 (49) donde el coeficiente de Manning [22] ´unicamente depende de la altura de rugosidad. 3.4. La aproximaci´on de Burguete En el modelo de Burguete [4], que considera flujo turbulento completamente desarrollado sobre la pared rugosa, el perfil de velocidades en el eje zse modela como: u(ζ) = (ulζ lbsi ζ > l 0 si ζ≤l(50) donde bes un par´ametro de ajuste que depende del ratio l/h y de las caracter´ısticas del flujo, la altura ζ=z−zbrepresenta la distancia vertical al fondo, lcaracteriza las rugosidades del fondo y la velocidad ules la velocidad a la distancia l. En este modelo, se asume que en canales lisos, la longitud l, puede identificarse el tama˜no de la subcapa viscosa, Figura 5. A partir de (50) es posible definir la velocidad promedio en la vertical: ¯u=1 hZζ=zb+h ζ=zb u dζ =ul b+ 1 h lb"1−l h1+b#(51) La misma hip´otesis de flujo turbulento en (50) nos permite escribir el esfuerzo en el fondo como τb=ρu2 l(52) 18
Figura 5: Perfil vertical de velocidad. donde es una constante adimensional dependiente de la rugosidad caracter´ıstica y del n´umero de Reynolds. Combinando (51) y (52) se obtiene la siguiente expresi´on para el coeficiente de fricci´on cf: cf=ρ(b+ 1)2l h2b"1−l h1+b#−2 (53) Este valor de coeficiente de fricci´on puede relacionarse con el obtenido en la aproximaci´on de Gauckler-Manning (48). Reorganizando (53) y escribi´endola frente a esta ´ultima: ρn21 h1/3≡(b+ 1)2l2b1 h2b"1−l h1+b#−2 (54) donde se ve que el perfil para la expresi´on tradicional de la fricci´on de Manning (48) es un caso concreto de (53) con b= 1/6, 1 h1/3≡1 h2b(55) con un nuevo t´ermino que modifica el valor de la fricci´on seg´un la relaci´on entre la distancia de rugosidad ly la altura h: "1−l h1+b#−2 (56) que est´a en acuerdo con leyes emp´ıricas para fondos rugosos que predicen un incremento para el coeficiente nde Manning para valores decrecientes de altura de l´amina de agua [17, 37]. Adem´as, se aprecia que gn2≡(b+ 1)2l2b(57) por lo que ahora n, a diferencia de la relaci´on (49), est´a caracterizada m´as en detalle por una distancia de rugosidad ly el exponente de velocidad b. En el perfil en (50) la longitud lse identifica el tama˜no de la subcapa viscosa. Considerando que la aproximaci´on de Manning-Strickler es un caso particular del perfil de velocidad m´as general en (39), asumiendo que K+ s≥70, la longitud 19
lpuede identificarse con el tama˜no de la subcapa rugosa en canales rugosos, [4]. En el caso de considerar Ks=l=d50, se puede aproximar como =gn2 (b+ 1)2l2b≈0,012 (58) Como se ve, el modelo de Burguete es capaz de explicar la aproximaci´on de Manning-Strickler interpretando el perfil de velocidad como una ley exponencial. Pero esta formulaci´on tiene puntos d´ebiles. El perfil de velocidad (50) no est´a formulado siguiendo los principios de similaridad completa en (30) y eso significa que necesita un exponente de ajuste b. En el modelo inicial en (42) el perfil sigue un patr´on logar´ıtmico en todos los casos, y ´unicamente aparece como inc´ognita la rugosidad equivalente Ks. El perfil (50) es discontinuo, ya que se asume la existencia de una subcapa de anchura ldonde la velocidad es despreciable y adem´as, cuando h→l, el esfuerzo τb,h→l=∞(59) es ilimitado. Adem´as, un perfil discontinuo de velocidades significa un perfil discontinuo de esfuerzos en la vertical, lo que contradice las observaciones experimentales. 20
4. Modelos de fricci´on en canales compuestos Es importante observar que muchos modelos emp´ıricos, como el modelo de Manning [22, 16], el modelo de Kellerhals [21] para rios de con lechos de gravas o el modelo de Darcy [7], fueron desarrollados para modelar p´erdidas de energ´ıa, asumiendo que el flujo se puede describir de forma adecuada a trav´es de una velocidad promedio ´unica en geometr´ıas unidimensionales. Las ecuaciones unidimensionales en flujo en canales de superficie libre o ecuaciones de Saint-Venant [31], pueden expresarse en forma conservativa como ∂U ∂t +∂F ∂x =H(60) donde U= (A, Q)TF= (Q, βQ2/A +gI1)T H= (0, g(I2+So)−Tb)T (61) Qes caudal circulante en la secci´on mojada A,So=−dz dx es la pendiente del fondo, y Tbes el esfuerzo cortante sobre toda la superficie de la secci´on. El coeficiente βaparece como resultante de asumir que existe una velocidad variable en la secci´on β=A Q2ZA u2dA (62) yI1eI2son las fuerzas de presi´on I1=ZH 0 σ(H−z)dz I2=ZH 0 ∂σ ∂x(H−z)dz (63) con σel ancho de la secci´on en la elevaci´on zyHel mayor valor de altura de l´amina de agua. En el modelo emp´ırico de Manning, desarrollado para flujo en canales unidimensionales, es com´un expresar la p´erdida de energ´ıa a trav´es de la pendiente de fricci´on adimensional Sf Sf=Tb gA =n2Q|Q| A2R4/3 h =n2U|U| R4/3 h (64) siendo Rhel radio hidr´aulico que se calcula como el ratio entre el ´area mojada Ay el per´ımetro mojado PyUla velocidad promedio en la secci´on. En el modelo bidimensional se puede escribir las p´erdidas de energ´ıa Sf,x ySf,y en los ejes xeyrespectivamente. Es destacable que la velocidad promedio Uen una secci´on compuesta por un cauce principal y una llanura de inundaci´on, se define como: U=Q A=Q Acauce +Allanura (65) aunque la los datos experimentales indican un gran variaci´on entre la velocidad del flujo en el cauce, Acauce y en la zona de la llanura Allanura. De esta forma el valor del coeficiente de Manning para ser correcto debe recoger la forma geom´etrica del canal, ya que la altura de la l´amina de agua cambia a lo largo de la secci´on [12, 11, 9]. 21
Cuadro 1: Valores del coeficiente cfpara los modelos MA,MB, y MC. Modelo cf MAgn2 h1/3 MBgn2 h1/3h1−l h7/6i−2 MCρ(b+ 1) l h2bh1−l h1+bi−2 Por lo tanto, es importante estudiar la aceptabilidad de las hip´otesis para el c´alculo del coeficiente de fricci´on n. Adem´as es importante comparar con modelos similares al modelo de Manning desarrollados utilizando una m´ınima cantidad adicional de informaci´on. Con ese fin, en esta secci´on se comparan tres modelos de fricci´on, que se llamar´an modelos MA,MB, y MC. El modelo MA es la aproximaci´on de Manning-Strickler en (48) y el modelo MBes un modelo intermedio que tiene en cuenta la longitud caracter´ıstica l, y no incorpora variaciones en el perfil de la velocidad, asumiendo b= 1/6 fijo e igual al modelo de Manning-Strickler. El modelo MCcontiene toda la informaci´on necesaria para desarrollar el perfil exponencial en (50). La Tabla 4 resume los valores del coeficiente cfpara cada uno de los modelos. Para ello se examinar´a la respuesta de los tres modelos frente a datos experimentales obtenidos en canales. 22
5. An´alisis de los datos experimentales Los experimentos Los datos experimentales utilizados en este trabajo provienen de medidas realizadas un canal modificable de 60 metros de largo por 10 metros de ancho (v´ease Fig. 6), en las instalaciones de la School of Civil Engineering de la Universidad de Birmingham [20]. En este canal se llevaron a cabo experimentos con distintas configuraciones geom´etricas con mayor o menor anchura de las llanuras adyacentes y la pendiente de separaci´on entre ellas. Figura 6: Fotos del canal y su esquema. De las distintas configuraciones posibles, se han elegido tres para este estudio, que se muestran en la Figura 7. 23
Figura 7: Geometr´ıa de los canales experimentales: (a) Canal I, (b) Canal II, (c) Canal III. Los datos se obtuvieron inyectando distintos caudales en cada canal, de los se han elegido cinco por cada canal como se detalla en la Tabla 2 marcados con una flecha, de tal manera que el calado Hcoincidiera entre los canales. Dicho calado Hes el que se proporciona en las bases de datos y es aqu´el que consegu´ıa la l´amina de agua en la secci´on en la que se midieron las velocidades en cada caso. Lo que no aparece en los datos, es el perfil de equilibrio de l´amina superficial, aunque se ha supuesto que se alcanzaba un estado estacionario. Para la nomenclatura de los casos se han respetado la de [20], se decir: n´umeros ascendentes seg´un el caudal, y algunos con guiones como segundas y terceras pruebas para ese mismo caudal: 7-1, 7-2..etc. Los datos Los datos medidos con los que se cuenta son: - La velocidad en direcci´on longitudinal, en funci´on de la elevaci´on y de la distancia al eje del cauce, esto es: ux(y, z) donde la posici´on xde medida se situ´o donde se daba un un cierto calado H que tambi´en se aporta. La Figura 8 muestra la representaci´on vectorial de las velocidades u(y, z) en la secci´on de altura H=0.199 m para mitad de canal I en el caso 5 (Q3=0.4511 m3/s). El perfil se considera sim´etrico en el eje yrespecto de la mitad del canal. - Una serie de par´ametros, algunos geom´etricos y conocidos ’a priori’ como la pendiente del 24
Cuadro 2: Casos con sus caudales y calados para los tres canales. Las flechas marcan los caudales Qjescogidos. Canal I Canal II Canal II Caso Q[m3/s]H[m]Caso Q[m3/s]H[m]Caso Q[m3/s]H[m] →1 0.20820 0.159 →1 0.21230 0.156 →1 0.22510 0.158 →2 0.23370 0.165 2-1 0.24920 0.170 →2 0.24120 0.166 3 0.28520 0.176 →2-2 0.24830 0.169 3 0.26760 0.177 4 0.35350 0.187 3 0.28210 0.178 4 0.30310 0.188 →5 0.45110 0.199 4-1 0.32382 0.187 →5 0.33230 0.198 6-1 0.60460 0.214 4-2 0.32370 0.187 →6 0.39220 0.215 →6-2 0.60010 0.214 →5 0.38320 0.198 7-1 0.55790 0.248 6-3 0.60010 0.214 →6 0.48000 0.214 →7-2 0.55810 0.248 7-1 1.01450 0.250 →7 0.76300 0.249 8-1 0.83600 0.299 →7-2 1.01450 0.250 8 1.11420 0.288 8-2 0.83490 0.300 7-3 1.01450 0.250 canal S0, las anchuras tanto del cauce central como de las llanuras de inundaci´on; otros medidos ’a posteriori’ como la pendiente de la superficie del agua Sw, estimaciones de n´umeros de Manning ´optimos y de n´umeros de Froude. El hecho de que los valores Sw≈S0indicaban que el perfil de equilibrio de la l´amina de agua era suave. Figura 8: Representaci´on vectorial de las velocidades ux(y, z) en la secci´on de altura Hpara mitad de canal I para Q3=0.4511 m3/s. Los datos de velocidad proporcionados ux(y, z) (Fig. 8) se pueden representar frente a la profundidad z, pero tambi´en promediar en la columna vertical de agua, dando lugar a una estimaci´on del promedio vertical ux(y), como se ver´a m´as adelante. 25
Figura 15: Adimensionalizaci´on de los perfiles de velocidad en el canal I, con Q=0.4511 m3/s . Figura 16: Parejas ´optimas (b,l) para los perfiles en las 48 columnas de agua para el caso 5 del canal I (v´ease Fig. 15). En rojo, los pares ´optimos para los perfiles en las llanuras, en negro los del cauce, en verde los de dos perfiles situados en la rampa que separa las dos zonas. 32
Figura 17: Error con los datos experimentales utilizando los valores adimensionales en (68) 33
6. Resultados mediante simulaci´on num´erica 2D La simulaci´on num´erica ofrece la oportunidad de comparar diferentes aspectos relacionados con la din´amica de los flujos en canales. Se va a analizar las tres aproximaciones detalladas en la Tabla 4, en el contexto de un modelo bidimensional, que permitir´a calcular en detalle las variaciones de velocidad en el eje perpendicular a la direcci´on principal del flujo. El esquema num´erico elegido es un m´etodo descentrado de primer orden en vol´umenes finitos. El esquema se detalla en la secci´on 9. Para obtener valores discretos de las variables en el espacio y en el tiempo es necesario dividir el dominio de integraci´on en celdas de c´alculo. Se ha elegido una malla rectangular estructurada caracterizada por ∆x= 0,3m y ∆y= 0,1m, con una longitud de 45 metros y una anchura que depende del canal. La altura del fondo de la malla depender´a de la forma del canal, seg´un se vio en los esquemas de la Fig. 7. Una representaci´on 3D de la variaci´on del fondo del canal I usando la malla descrita se ve en la Fig. 18. 0 0.05 0.1 0.15 0.2 Z 0 10 20 30 40 X 0246810 Y XY Z Z 0.1838 0.1716 0.1593 0.1471 0.1348 0.1226 0.1104 0.0981 0.0859 0.0736 0.0614 0.0491 0.0369 0.0246 0.0124 Figura 18: Representaci´on 3D de la variaci´on del fondo del canal I, utilizando una malla rectangular: ∆x= 0,3m y ∆y= 0,1m. La elecci´on del tipo de mallado viene dada por la necesidad de caracterizar adecuadamente las variaciones espaciales de las variables de c´alculo que son mucho mayores en la direcci´on transversal comparadas con las variaciones en la direcci´on longitudinal. Una malla de densidad constante no ofrece en situaciones como las estudiadas una mejor resoluci´on de las variables de c´alculo pero s´ı que exige un coste computacional m´as alto. La Figura 19 muestra un detalle tridimensional de la variaci´on del fondo del canal II, utilizando una malla rectangular: ∆x= 0,3m y ∆y= 0,1m. Condiciones de estacionariedad El sistema debe llegar al estado estacionario: s´olo ah´ı las variables tendr´an sentido, pues se quieren emular las condiciones en las que se midieron los datos experimentales en el canal de laboratorio. En cada uno de los experimentos se midieron las velocidades en una zona del canal con un calado Hdeterminado y por tanto, ahora en la simulaci´on tambi´en se ha de encontrar esa secci´on 34
0 0.1 0.2 0.3 0.4 Z 40 0246810 Y XY Z Z 0.4182 0.3903 0.3625 0.3346 0.3067 0.2788 0.2510 0.2231 0.1952 0.1674 0.1395 0.1116 0.0838 0.0559 0.0280 Figura 19: Detalle tridimensional de la variaci´on del fondo del canal II, utilizando una malla rectangular: ∆x= 0,3m y ∆y= 0,1m. transversal que tenga precisamente ese calado. Una vez encontrada esa secci´on, se comparar´an las velocidades experimentales con las calculadas seg´un el esquema num´erico. Condiciones de contorno Utilizar un modelo de fricci´on u otro implica caracterizar la celda con un valor de cfdistinto. Esto conlleva diferencias en la altura de la l´amina de agua, diferencias energ´eticas y distintas velocidades. Se podr´ıa pensar que no se producen pr´acticamente cambios en el calado a lo largo del canal y que las velocidades tampoco cambian sustancialmente. En ese caso bastar´ıa forzar un calado en la salida igual al Hque se busca, y ya servir´ıan los datos de velocidad en cualquier secci´on. Y la condici´on de contorno en la entrada es el caudal Qinyectado en el experimento. Se han calculado en simulaci´on la variaci´on del calado para distintos valores del n´umero de Manning con el objetivo de tener una idea de c´omo debe ser la condici´on de contorno en la salida del canal, pues convendr´ıa que se acercase el nivel de la l´amina de agua al valor de H. Los resultados para los tres canales elegiendo su caudal Q3respectivo se pueden ver en la Fig. 20, donde se ha elegido un hipot´etico canal de 300 metros aunque en lo sucesivo se vaya a prestar atenci´on a 35
los 45 ´ultimos metros. Se intuye que con el valor del n´umero de Manning ´optimo recomendado por [20], n = 0.00904 para este caso, el perfil quedar´a por debajo de aqu´el con n=0.01000 y habr´a que escoger un calado en la salida mayor que el buscado H=0.199 m para llegar a ´el. Para este caso se elegi´o un valor de 0.210 m. Si se representan las velocidades promediadas en la vertical sobre la columna de agua central en los tres canales, se llega a la Fig. 21. Se puede apreciar que en las posiciones donde las curvas cortan con H= 0,199 m la velocidad no se acerca a la experimental en el centro del canal (ucentro=0.86 m/s). Sin embargo, se sabe que el flujo es lento, subcr´ıtico seg´un los valores de Froude ofrecidos (Fr< 1), por lo que no se deber´ıa medir donde el perfil no sea suave, pues en esas zonas no se estar´a modelando el flujo fielmente. V´ease en la Fig. ?? que el perfil de la l´amina libre no parece constante en la longitud que ocupar´ıa el canal (45 m). Otra opci´on es no forzar ning´un calado concreto a la salida, permitiendo que el propio flujo sea el que lo establezca. De esta forma se ha conseguido un perfil de equilibrio constante y es la forma en la que se calcular´a en lo sucesivo. Fijando una salida o no, en todo caso lo que se ha buscado es llegar a un calado constante, que es el s´ıntoma de haber llegado al estado estacionario. 36
Figura 20: Variaci´on del calado para distintos n´umeros de Manning con la salida forzada a h = 0.210 m 37
Figura 21: Variaci´on de la velocidad promediada sobre el centro del canal para distintos n´umeros de Manning, con el calado a la salida forzado a h = 0.210 m 38
B´usqueda de la secci´on de medida Como ya ha quedado patente, en la variaci´on del calado a lo largo del canal es decisiva la fricci´on, de tal manera que ese calado Hdonde fueron medidos los datos, no ocupa siempre el mismo lugar en la longitud del canal. Un vez encontrada esa secci´on, habr´an de compararse las velocidades experimentales con las calculadas por el esquema num´erico, en concreto los promedios en la vertical sobre la columna de agua. En la Fig. 22 se muestran las velocidades experimentales promedio para los distintos caudales contra los que habr´a que comparar los resultados num´ericos. Como el perfil de equilibrio tiende a ser constante en una zona, realmente es l´ıcito medir la velocidad en cualquier secci´on de esa zona, m´as que buscar una secci´on concreta. De hecho se hace referencia a una zona de test en los datos de [20], y no a una secci´on determinada. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Velocidad u(y) promediada en z [m/s] Distancia a mitad de canal [m] Figura 22: Velocidades experimentales promedio para los distintos caudales en la Tabla 5 en el canal I. Q5=1.0145 m3/s (4), Q4=0.6001 m3/s (•), Q3=0.4511 m3/s (◦), Q2=0.2082 m3/s (). Volviendo a los resultados de simulaci´on del caso del canal I y caudal Q3=0.4511 m3/s, se han le´ıdo los valores de velocidad en ciertas partes del canal definidas de antemano (”sondas”) y que se han dispuesto equiespaciadas: a 0, 9, 18, 27, 36 y 45 metros de la entrada. De esta forma, para cada simulaci´on con un modelo distinto, se ”leer´a”la sonda que m´as cerca se encuentre del lugar donde los calados num´ericos se cortan con el calado experimental H=0.199 m y se comparar´a con la curva experimental. 39
Resultados As´ı, se han simulado los tres modelos propuestos, para el caso 5 (Q3= 0,4511m3/s) del canal I, con los siguientes par´ametros que se muestran en la Tabla 6 Cuadro 4: Simulaciones en el canal I para el caso 5 y los tres modelos. Experimento num´erico Modelo usado Par´ametros 1MAn= 0,00904 2MBn= 0,00904, l = 0,00001 3MBn= 0,00904, l = 0,00010 4MBn= 0,00904, l = 0,00100 5MBn= 0,00904, l = 0,01000 6MCe= 0,01250, b = 0,15, l = 0,00100 7MCe= 0,01250, b = 0,15, l = 0,01000 8MCe= 0,01250, b = 0,14, l = 7,5∆10−5 donde la nelegida es la que suger´ıa [20], , fue calculada en (??) y se ha optado por utilizar distintos valores de l, de tal forma que uno fuera diez veces mayor que el anterior. Para el experimento 8 se eligieron bylcomo el par ´optimo que minimiza el error total, como se vio en el an´alisis de los modelos y en concreto en la Fig. 17. Debe asegurarse que la l´amina de agua llegue a una altura razonablemente constante y observar de qu´e calado se trata (Fig. ??). Se advierte que cuanto m´as se ha acercado el modelo al calado Hbuscado (H=0.199 m en el caso de Q3=0.4511 m3/s), m´as se asemeja la velocidad en el centro del canal a la experimental (v´ease Fig. ??). Si se quieren comparar las velocidades a lo largo de la secci´on, entonces hay que fijarse en la Fig. 25 donde se aprecia que el modelo MCes el que mejor se ajusta al caso, y m´as en detalle en la Fig. 26 donde se muestran las velocidades en la ”sonda.a9 metros de la entrada. De esta manera, para el modelo MBse sigue manteniendo este n, pero ahora se le ha a˜nadido una rugosidad l. En realidad, los modelos MAyMBpara rugosidades (l) peque˜nas se comportan pr´acticamente igual, y para valores grandes ya empeora. Esto estar´ıa de acuerdo con que las rugosidades deber´ıan ser mucho m´as peque˜nas que el calado (l << h), por lo que valores de lde cent´ımetros ser´ıan demasiado grandes para una situaci´on donde en la llanura se van a tener calados de alrededor de 5 cm, pues ya no se cumplir´ıa la hip´otesis anterior. Si se desea ver qu´e tal predicen los tres modelos, para un mismo canal se puede variar el caudal y observar las nuevas velocidades y calados, y compararlos con los experimentales. Si se recuerda la Tabla 4, el valor del coeficiente cfva teniendo en cuenta m´as par´ametros a la hora de calcular, si bien los tres modelos tienen en cuenta la altura de la l´amina de agua. Por ejemplo, se ha inyectado ahora en el canal I un caudal Q5= 1,0145m3/s, pero se siguen usando los par´ametros ´optimos para el caso con Q3= 0,45110m3/s. Los resultados se se pueden observar en las Figs. 27-30. Se aprecia que los mismos par´ametros que funcionan bien para un caso, 40
Figura 23: Calado a lo largo del canal I para el caso 5 (Q3=0.4511 m3/s) seg´un el c´alculo con distintos modelos no tienen por qu´e hacerlo cuando se cambian las variables hidrodin´amicas del problema (caudal y calado). El modelo MCten´ıa en cuenta argumentos energ´eticos, como se resum´ıa en la expresi´on ??, y su comportamiento ya no tendr´a por qu´e parecerse a MAni a MB. Esto queda patente en la figura anterior ??, donde se aprecia lo cercanos (superpuestos incluso) que aparecen los calados para los modelos MAyMBpero lo separados que se encuentran de MC. En el modelo C intervienen una serie de factores (bo exponente del perfil, o coeficiente aerodin´amico que es un t´ermino de energ´ıa, y lo rugosidad equivalente del fondo) que se relacionan de manera compleja y la relaci´on con altura del calado, por ejemplo, ya no ser´a predecible como lo era con el par´ametro nde Manning en la Fig. 20. Pese a la mayor dificultad de este modelo, en este trabajo ha sido el que mejor ha ajustado la velocidad calculada num´ericamente a la experimental. 41
[16] P. G. Gauckler.´ Etudes th´eoriques et pratiques sur l’´ecoulement et le mouvement des eaux. Comptes Rendues de l’Acad´emie des Sciences, Paris.1867. [17] R. D Jarrett. Hydraulics of high-gradient streams. J. Hydraul.Eng.,110, 1519–1539, 1984. [18] R.J. Leveque.Finite Volume Methods fot Hyperbolic Problems. Cambridge University Press, New York, 2002, p. 311, 2002. [19] R.M. Ferreira.River morphodynamics and sediment transport conceptual model and solutions. Phd. Thesis, Lisbon, 2005. [20] D. W. Knight. FCF Series. Flow Database at the Univ. of Birmingham. www.flowdata.bham.ac.uk 2007. [21] R. Kellerhals. Stable channels with gravel-paved beds. J. Wtrwy. and Harb. Div.,93, 63–84, 1967. [22] R. Manning. On the flow of water in open channels and pipes. Trans. Inst. Civil Engineers, 20, 161–207,1890. [23] A. S. Monin and A. M. Yaglom (1975) Statistical fluid mechanics. Vol. 2: Mechanics of turbulence. MIT Press,1975. [24] J. Murillo, J. Burguete, P. Brufau, and P. Garc´ ıa-Navarro. The influence of source terms on stability, accuracy and conservation in two-dimensional shallow flow simulation using triangular finite volumes. International Journal of Numerical Methods in Fluids 54, 543–590, 2007. [25] J. Murillo, P. Garc´ ıa-Navarro and J. Burguete. Time Step Restrictions For Well Balanced Shallow Water Solutions In Non-Zero Velocity Steady States. International Journal of Numerical Methods in Fluids 56, 661–686, 2008. [26] J. Murillo, P. Garc´ ıa-Navarro and J. Burguete. Conservative Numerical Simulation of Multicomponent Transport in Two-Dimensional Unsteady Shallow Water Flow, Journal of Computational Physics,228, 5539–5573,2009. [27] J. Murillo, P. Garcia-Navarro. Weak solutions for partial differential equations with source terms: Application to the shallow water equations,Journal of Computational Physics 229, 4327–4368, 2010. [28] M. R. Raupach, R. A. Antonia, and S. Rajagopalan .Rough-wall turbulent boundary layers, Appl. Mech. Review,44, 1–25,1991. [29] L.C. van Rijn.Sediment Transport, Part III: Bed forms and Alluvial Roughness, J. Hydr. Engrg.,110, 1733–1740, 1984. [30] P.L. Roe.A basis for Upwind Differencing of the Two-Dimensional Unsteady Euler Equations. Numerical Methods in Fluid Dynamics, Vol II. Oxford University Press, Oxford, 1986. 48
[31] A. J. C. B. de Saint-Venant. Th´eorie de mouvement nonpermanent des eaux avec application aux crues de rivi`eres et ´a l’introduction des mar´ees dans leur lit. Comptes Rendues de l’Acad´emie des SciencesComptes Rendues de l’Acad´emie des Sciences. Paris, 1871. [32] E.F. Toro.Shock-Capturing Methods for Free-Surface Shallow Flows. Wiley, New York, 2001, p. 109, 2001. [33] A. A. Townsend.The Structure of Turbulent Shear Flows.2nd edition, Cambridge Univ. Press.1976 [34] M.E. V´ azquez-Cend´ on. Improved treatment of source terms in upwind schemes for the shallow water equations in channels with irregular geometry. Journal of Computational Physics 148, 497–498, 1999. [35] C.B. Vreugdenhil.Numerical methods for shallow-water flow. Kluwer Ac. Pub., Dordrecht, The Netherlands, 1994. [36] F. M. White.Fluid Mechanics.. McGaw-Hill International Editions,New York, 1986. [37] F. M. White.Viscous fluid flow. McGaw-Hill International Editions,New York, 1991. [38] W. Wu. Depth-Averaged Two-Dimensional Numerical Modeling of Unsteady Flow and Nonuniform Sediment Transport in Open Channels. J. Hydr. Engrg. 130, 1013–1020, 2004. 49
. 50
9. Anexo: El esquema num´erico Las ecuaciones de aguas poco profundas se deducen a partir de la ecuaci´on de conservaci´on de la masa y de las ecuaciones de Navier-Stokes para flujo incompresible. Supuestas la velocidad en la escala vertical, w, mucho menor que en las direcciones del plano x, y y la presi´on, hidro´estatica, el proceso de promediar en la vertical las ecuaciones convierte el problema tridimensional en uno bidimensional de grosor variable h[15]. El conjunto de ecuaciones que describen lan conservaci´on de la masa y la conservaci´on del momento generan dan lugar a un sistema de ecuaciones en derivadas parciales se puede expresar de forma compacta como sigue [26]: ∂U ∂t +∂F(U) ∂x +∂G(U) ∂y =S(U) (69) donde U=h qxqyT(70) son las variables conservadas, con hque representa la profundidad de la columna de agua y qxy qylas componentes del vector velocidad, ua lo largo del eje x e y respectivamente y promediadas en la vertical. Los flujos de estas variables vienen dadas por: F=qx q2 y h+1 2gh2qxqy/h T G=qyqxqy/h q2 y h+1 2gh2T(71) donde ges la aceleraci´on de la gravedad. Los t´erminos fuente del sistema se dividen en tres tipos. Los t´erminos que hacen referencia a la pendiente del fondo y a la fricci´on en las ecuaciones de momento son: S=0gh(Sox −Sfx)gh(Soy −Sfy)T(72) donde Sox =−∂z ∂x Soy =−∂z ∂y (73) son la pendientes del fondo en los ejes xey, definidas a trav´es del nivel zySfx ySfy representan las p´erdidas por fricci´on. El sistema (69) depende del tiempo y no lineal. Bajo la hip´otesis de dominio convectivo, el sistema se puede clasificar y tratar desde el punto de vista num´erico como perteneciente a la familia de los sistemas hiperb´olicos. Las propiedades matem´aticas de (69) incluyen la existencia de una matriz Jacobiana Jn, del flujo normal En =Fnx+Gny, definida como: Jn=∂En ∂U=∂F ∂Unx+∂G ∂Uny(74) Esta matriz Jacobiana ser´a la base para llevar a cabo la discretizaci´on num´erica que se presenta en este trabajo. Para comenzar la t´ecnica de vol´umenes finitos, el sistema (69) se integra en el volumen o malla de celdas, Ω : 51
Figura 31: Representaci´on constante de las varaibles en cada celda. ∂ ∂t ZΩ UdΩ + ZΩ (−→ ∇E)dΩ = ZΩ SdΩ (75) donde se asume que la tercera integral se puede formular de la siguiente manera [25, 26, 27]: ZΩ SdΩ = I∂Ω (Tn)dl (76) donde Tes una matriz adecuada. Esto lleva a la siguiente formulaci´on: ∂ ∂t ZΩ UdΩ + I∂Ω Endl =I∂Ω Tndl (77) Cuando el dominio se subdivide en celdas Ωien una malla fija en el tiempo, Figura 31, la ecuaci´on (77) tambi´en se puede aplicar a cada celda. ∂ ∂t ZΩ UidΩi+ NE X k=1 Zek+1 ek Ejnklk= NE X k=1 Zek+1 ek Tnklk(78) En primer orden las cantidades vectoriales son uniformes en cada celda y la ecuaci´on (77) queda reducida a: (Un+1 i−Un i) ∆tAi+ NE X k=1 (δE−T)knklk= 0 (79) donde δE=Ej−Ei, con EjyEiel valor de la funci´on Een la celda vecina jy en la celda i respectivamente y conectadas a trav´es del lado k,nkes el vector normal hacia afuera del borde de la celda k,lkes la longitud correspondiente a dicho borde, NE es el n´umero de bordes que necesita para definir la celda y Tk Debido al car´acter no lineal del flujo E, la aproximaci´on a un flujo Jacobiano, e Jn,k [30] permite una linealizaci´on local dando lugar a un sistema de 3 + pvalores propios reales e λm ky vectores propios eem k, que se construyen con las siguientes variables promedio [24, 26, 27] euk=ui√hi+uj√hj √hi+√hjevk=vi√hi+vj√hj √hi+√hjeck=qghi+hj 2(80) 52
dando lugar a e λ1 k= (e un +ec)ke λ2 k= (e un)ke λ3 k= (e un −ec)k(81) y ee1 k= 1 eu+ecnx ev+ecny kee2 k= 1 −ecny ecnx kee3 k= 1 eu−ecnx ev−ecny k (82) Las matrices e Pk, t e P−1 k, se pueden construir a partir de los vectores propios eem kde e Jn,k de forma que la diagonalicen e Jn,k = (e PΛe P−1)ke Pk=ee1 kee2 kee3 k(83) donde e λm kson los valores propios de la matriz diagonal Λ. Tambi´en, a partir de la matriz Jacobiana aproximada [30] e Jn,keem k=e λm keem km= 1,2,3,...,3 + p(84) El problema se reduce a un problema unidimensional proyectado sobre la direcci´on nen cada borde de celda [32]. Adem´as, la diferencia del vector Ua trav´es de cada lado de la celda se proyecta sobre la base de vectores propios δUk= 3 X m=1 (αee)m k(85) Donde las expresiones que dan los coeficientes αkson: α1,3 k=δhk 2±1 2eck(δqk−e ukδhk)nkα2 k=1 2eck(δqk−e ukδhk)nT,k (86) con nT,k = (−ny, nx). La contribuci´on de δ(En)ken una celda kse puede escribir como: δ(En)k= 3 X m=1 (e λαee)m klk(87) La pendiente del fondo y el t´ermino de fricci´on se pueden descomponer en la base de los vectores propios con el objetivo de asegurar el equilibrio discreto con los t´erminos de flujo (87), de tal manera que se asegura en los casos estacionarios con velocidad nula y no nula [27, 34]: (Tn)k=e PkBklk= 3 X m=1 (βmeem)klk(88) con Bk=β1β2β3T k. Los coeficientes son β1,3 k=∓ eck 2(δz +dnSf)kβ2 k= 0 (89) El esquema descentrado expl´ıcito de primer para un sistema no-reactivo y no-difusivo toma la forma 53
Figura 32: Selecci´on de la informaci´on necesaria en el m´etodo descentrado. Un+1 i=Un i+ ∆t NE X k=1 Ψn i,k (90) donde la contribuci´on de cada lado de la celda, Ψi,k, s´olo recoge la informaci´on en la direcci´on entrante, Figura (32): Ψi,k = 3 X m=1 ((e λ−α−β−)ee)m klk/Ai(91) donde e λ−=1 2(e λ−|e λ|) y β−=1 2(β−|β|). 54