Este documento presenta una nueva familia de esquemas híbridos semi-implícitos de volumen finito/elemento finito en mallas escalonadas basadas en bordes para la solución numérica de las ecuaciones de Reynolds-Averaged Navier-Stokes (RANS) incompresibles en combinación con el modelo de turbulencia. La reología para el cálculo del coeficiente de viscosidad laminar considerado en este trabajo es la de un fluido Herschel-Bulkley no newtoniano (ley de potencia) con esfuerzo de cedencia, que incluye al fluido de Bingham y a los fluidos newtonianos clásicos como casos especiales. Para la discretización espacial, utilizamos mallas simples no estructuradas escalonadas basadas en bordes, así como mallas cartesiana no uniformes escalonadas. Con el fin de obtener un algoritmo simple y computacionalmente eficiente, aplicamos una técnica de división de operadores, donde los términos convectivos hiperbólicos de las ecuaciones RANS se discretizan explícitamente mediante un esquema de volumen finito tipo Godunov, mientras que los términos parabólicos viscosos, los términos de presión elípticos y los términos de fuente algebraicos rígidos del modelo se discretizan de forma implícita. Para la discretización de la ecuación de Poisson de presión elíptica, utilizamos elementos finitos clásicos y conformes en triángulos y rectángulos, respectivamente. La discretización implícita de los términos viscosos es obligatoria para los fluidos no newtonianos, ya que la viscosidad aparente puede tender a infinito para los fluidos con esfuerzo de cedencia y ciertos fluidos de ley de potencia. Se lleva a cabo con elementos finitos en mallas simples triangulares y con volúmenes finitos en rectángulos. Para mallas cartesianas y mallas no estructuradas ortogonales más generales, podemos demostrar que nuestro nuevo esquema puede preservar el de y . Esto se logra mediante una discretización implícita especial de los términos de fuente de relajación algebraica rígida, utilizando una combinación adecuada de las ecuaciones de evolución discretas para los logaritmos de y . El método se aplica a algunos problemas de referencia académica clásicos para flujos no newtonianos y turbulentos en dos dimensiones espaciales, comparando los resultados numéricos obtenidos con soluciones de referencia exactas o numéricas disponibles. En todos los casos, se observa una excelente concordancia.