El modelo simula la interacción de tres procesos: el comportamiento del viento, la trayectoria de las partículas en salto y la transferencia de momento entre el flujo de viento y las partículas suspendidas. Es un modelo en tres dimensiones que incluye flujo turbulento k-e, fuerza arrastre, gravedad e interacciones entre partículas y fluido. El modelo rechaza la colisión entre partículas, las fuerzas electrostáticas y el giro de partículas sobre su propio eje (Fuerza de Magnus) ya que se busca realizar un modelo simplificado, que demuestre un comportamiento idealizado de saltación. Sin embargo, se podría incorporar este tipo de fuerzas incluyendo una física que gobierne estos efectos. El modelo turbulento utilizado es una adaptación del modelo estándar de k-e tomando en cuenta constantes sugeridas por Richard & Hoxey [8].

COMSOL MULTIPHYSICS 5.1® es utilizado para la creación de geometrías, generación de mallas, simulación el flujo turbulento de aire y el salto de partículas. Las simulaciones se realizaron en un sistema operativo Windows con procesador Intel® Core™ i5-4200U CPU @ 160 GHz.

El dominio de estudio es una caja que se extiende 2 m en el eje X, 0.5 m en Y y 0.5 m en Z como se muestra en la Figura 2. Adicionalmente, se generaron 10 capas de 0.01 m de alto desde la base de la caja con el fin de obtener una mejor aproximación de flujo turbulento y del fenómeno de salto cerca del terreno (base) de la caja. La geometría utiliza un tipo de mallado de tetraedros libres ( free tetrahedral ) como se muestra en la Figura 3. En las capas superiores se estableció una malla normal ( normal ) y fina ( fine ), como se muestra en la Figura 3a y 3b, respectivamente. En las 5 capas inferiores se estableció una malla extra fina (extra fine) por ser una zona crítica donde se genera la mayor cantidad de saltos (Figura 3c). Comsol es un programa iterativo que según el mallado genera diversos elementos para mejorar la solución final. Este tipo de mallado generó un total de 680 495 elementos en todo el dominio.

Geometría del dominio de estudio. La caja se extiende 2 m en X, 0.5 m en Y y 0.5 m en Z.

El perfil de viento es clave en el salto de partículas pues regula la cantidad de momento que se transfiere a la partícula en salto. La ecuación de perfil de viento adoptada para las simulaciones es logarítmica y en estado estacionario representada por la ecuación 1 [9-10].

(1)

Donde z [m] es la distancia vertical desde la base, u* [m/s] es la velocidad de cizalladura del viento, k es la constante de Von Karman y z_o [m] es la rugosidad de la superficie. La velocidad de cizalladura del viento puede ser calculada según la ecuación 2 [9-10].

(2)

Donde U_h (m/s) es la velocidad de referencia a la altura de h [m].

Para el desarrollo del modelo en Comsol se utilizó el módulo de transporte de fluidos con la física de Turbulent Flow, k.e (spf). Este es un modelo de dos ecuaciones que sirve para simular las características de un fluido turbulento. La primera variable de transporte, k, determina la energía cinética turbulenta y la segunda variable, e, determina la disipación de energía turbulenta. Ambas variables están determinadas por las ecuaciones 3 y 4 [8].

(3)

Donde Cu es una constante del modelo igual a 0.013 [8].

(4)

Para la simulación de partículas en salto se utiliza una velocidad de cizalladura del viento de 0.0441 m/s. Esta velocidad, a una altura de 10 m, es equivalente a 1 m/s y según la escala de Beautfort se describe como una ventolina en la cual el humo indica la dirección del viento [11]. Se utiliza la constante de Von Karman con un valor de 0.42 y un valor de de 7.31x10^-4 m que representa un valor equivalente a la rugosidad aerodinámica típica deun desierto [7].

Para la generación del perfil de viento se realizó un estudio en estado estacionario donde se definió una cara de ingreso ( inlet ) con la condición de velocidad en el eje X equivalente a la ecuación 2; una cara de salida ( pressure oulet ), con una condición de borde de presión igual a 1 atm; la cara superior se estableció una cara de frontera abierta donde no exista estrés viscoso; y, en las caras laterales se impuso una función de simetría, la cual describe que lo que sucede al interior de las paredes laterales sucede a su exterior. Dentro de los ajustes de terreno se utiliza una función de pared ( wall functions ), en la cual se aplica una rugosidad equivalente, , y este valor se calcula con la ecuación 5 [12].

(5)

El movimiento de las partículas en salto está determinado principalmente por la fuerza gravitacional y las fuerzas proporcionadas por el fluido propiamente. La principal fuerza que ejerce el fluido y afecta la trayectoria es la fuerza de arrastre ( drag force ), F_d . Esta fuerza está determinada por las ecuaciones 6, 7 y 8 [13].

(6)

(7)

(8)

Donde τ_p [s] es el tiempo de respuesta de la velocidad de partícula, m _p [kg] es la masa, v [m/s]es la velocidad de las partículas, es la viscosidad dinámica del fluido [Pa·s], ρ_p [kg/m3] es la densidad de la partícula, d_p [m] diámetro de la partícula y u’ [m/s] es la velocidad del fluido. Al utilizar un flujo con disipación turbulenta, el campo de velocidad del fluido se adapta a la ecuación 8.

Esta ecuación consta de dos términos: primero u [m/s] representa la velocidad media del fluido y el segundo término representa una fluctuación turbulenta donde k [m2/s2] es la energía cinética turbulenta y ζ es un número aleatorio distribuido normalmente con media de cero y desviación estándar de uno generado por el programa. Estas ecuaciones son adaptaciones de la ley de arrastre establecida por George Gabriel Stokes con el modelo de flujos cargados de partículas. Este modelo se refiere a una clase de flujo de fluido de dos fases en el que una de las fases está conectada de forma continua (fase continua o vehículo) y la otra fase se compone de partículas pequeñas inmiscibles (fase dispersa o de partículas) [14].

Las partículas saltantes al encontrarse en el sistema gravitatorio terrestre “sienten” la fuerza de gravedad Fg determinada por la ecuación 9 [13].

(9)

Donde g [m/s2] es la aceleración de la gravedad y ρ [kg/m3] es la densidad del aire.

Para determinar las trayectorias se realizó un estudio dependiente del tiempo en intervalos de 0.01-0.001 s, dependiendo del estudio a realizarse. Por ejemplo, para obtener una mejorvisualización de los saltos realizados a lo largo de la superficie, se debe asignar un menor intervalo de tiempo para poder captar el momento exacto en que la partícula topa el suelo. En las propiedades de las partículas se fijó una densidad promedio de polvo igual a 2.5 g/cm³ y se realizó el estudio para diámetros de partículas de 100, 250 y 500 mm. Se ingresaron 20 partículas dispersas por toda la cara frontal a diferentes alturas que serán arrastradas de acuerdo a las ecuaciones que representan la trayectoria de partículas y el perfil de viento. En el terreno se mantiene la rugosidad establecida por el perfil de viento y se adiciona una opción de rebote. Las condiciones para la fuerza de arrastre como el campo de velocidad y la viscosidad dinámica van a depender del estudio estacionario de perfil de viento.

a) Corte de velocidades en el dominio de estudio. b) Corte A1 para determinación del perfil de viento.

Para el estudio de salto de partículas se arrojaron 20 partículas a diferentes alturas sin velocidad inicial, esto quiere decir que ganarán momento netamente del perfil de viento descrito previamente en los métodos y materiales. En la Figura 6a) se aprecia una representación de las 20 partículas en su posición inicial en la cara frontal. Las partículas están coloreadas en función a la velocidad que el viento transferirá de acuerdo a la altura.

Ya que el terreno es casi ideal, no hay geometrías que obstaculicen la trayectoria y la velocidad del viento ingresa linealmente a lo largo del eje X. Las partículas de estudio casi no“sentirán” desviaciones en el eje Y. No obstante, este modelo puede ser aplicado en dominios que sí incluyan efectos de obstrucciones y topográficos.

En la Figura 6b-d se observa la posición de las partículas después de 3 s, en este punto se nota que todas las partículas están casi a nivel del suelo y con la mayor parte de su energía disipada.

a) Posición inicial de partículas, b) Posición después de 3 (s), partículas de 100 mm de diámetro, c) Posición después de 3 (s), partículas de 250 mm de diámetro, d) Posición después de 3 (s), :partículas de 500 mm de diámetro.

El modelo propuesto se simuló tomando las condiciones y ajustes de contornos, tanto para el perfil de viento como para el salto de partículas, establecidas por Jasper F. Kok y Nilton O. Renno descritas en el modelo COMSALT [7]. El mencionado estudio utilizó una velocidad de cizallamiento igual a 0.4 m/s y una rugosidad equivalente representada por la ecuación 10 [7,10,15].

Además, las partículas fueron expulsadas desde la superficie con una rapidez inicial de 1 m/s y un ángulo de 40o desde la horizontal. Los datos obtenidos y la comparación con respecto al modelo COMSALT se representan en la Figura 8.

(10)

Las trayectorias obtenidas en el modelo descrito en este trabajo de investigación muestran un comportamiento semejante a las trayectorias obtenidas en el modelo COMSALT [7] que fueron validadas por una serie de datos experimentales. Las partículas más pequeñas descritas en la Figura 8a) tienden a viajar distancias mayores en el eje horizontal. Por otro lado, las partículas de mayor tamaño (Figure 8b) ganan más momento vertical, pero al tener una mayor masa son mayormente afectadas por la fuerza de gravedad por lo que tienden a caer más rápido.

Como se mencionó anteriormente, las partículas con mayor diámetro son las causantes que partículas menores entren en salto. En la Figura 7a) se muestra cómo estas partículas generan la mayor cantidad de saltos a lo largo de una superficie, ya que, al obtener una mayor cantidad de momento vertical, generan una mayor energía potencial lo que conlleva a que puedan realizar una mayor cantidad de saltos a lo largo de la superficie en comparación a partículas con menor masa. Además, se realizó un análisis de conservación de energías en diferentes puntos para cuantificar la energía perdida en cada rebote. La ecuación para determinar la energía perdida está dada por la ecuación 11[16].

(11)

Donde E_Ma y E_Mb es la energía mecánica en el punto a y b respectivamente, los puntos a y b son elegidos aleatoriamente entre cada rebote que realiza la partícula, y W_FNC es el trabajo de las fuerzas no conservativas que en este caso será la energía perdida por fricción en cada contacto de la partícula con la superficie. La energía mecánica está representada por la suma entre la energía cinética ( E_k ), ecuación 12, y potencial ( E_p ),ecuación 13, en cada punto respectivamente [16].

(12)

(13)

Donde v es la velocidad de la partícula y h la altura en la que se encuentra la partícula en un momento determinado.

Se realizó este estudio de conservación de energía tomando en cuenta el antes y después de cada salto para la partícula 19, representada en la Figura 7b), ya que es la que más energía potencial tiene y da la mayor cantidad de saltos a lo largo de la superficie. En cada salto realizado por la partícula hay una pérdida promedio del 15± 5% en cada salto.

Comparación de trayectorias de partículas con diámetros de (a) 100 μm, (b) 250 μm. En el lado derecho se presentan trayectorias obtenidas con el modelo de Comsol y en el lado izquierdo los resultados por COMSALT. Las figuras se han reproducido con los permisos de autor Kok & Renno (2009) [4].

Con el fin de mostrar las bondades del modelo, en la Figura 9 se muestra un dominio con terreno irregular y con cortes de velocidad del viento y la trayectoria de una partícula de 500 μm de diámetro. Como es evidente, la topografía juega un papel importante en el salto de las partículas desviándolas de su comportamiento ideal, es decir, cuando el terreno es plano.

En la Figura 9b) se nota claramente que el viento pierde su tendencia logarítmica. En las zonas con mayor irregularidad existe una aceleración del fluido para poder superar el obstáculo generado por el propio terreno. Esto causa que las partículas que pasen por estas zonas, experimenten una mayor aceleración, lo cual causara desviaciones en la trayectoria de estas. Por otro lado, en la Figura 9c) muestra la trayectoria de una partícula de 500 μm de diámetro. Esta partícula es transportada a lo largo de la superficie por una serie de saltos, con movimientos hacia adelante y hacia atrás. Este tipo de movimiento ocurre por la irregularidad del terreno, el tipo de salto va a depender directamente del ángulo de impacto sobre la superficie inclinada [17].

(a) Geometría con terreno irregular, (b) corte de velocidades, (c) trayectoria de partícula de 500 μm de diámetro.