Método Falsa Posición: Calcular factor de fricción

PROBLEMA: En una sección de tubo, la caída de presión se calcula así:

$$\Delta p=f\frac { L\rho { V }^{ 2 } }{ 2D }$$

donde Dp = caída de presión (Pa), f = factor de fricción, L = longitud del tubo [m], r = densidad (kg/m³),

V = velocidad (m/s), y D = diámetro (m). Para el flujo turbulento, la ecuación de Colebrook proporciona un medio para calcular el factor de fricción,

$$\frac { 1 }{ \sqrt { f }  } =-2{ log }_{ 10 }\left( \frac { \varepsilon  }{ 3.7D } +\frac { 2.51 }{ Re\sqrt { f }  }  \right) $$

donde e = rugosidad (mm), y Re = número de Reynolds,

$$Re=\frac { \rho VD }{ \mu  }$$

donde m = viscosidad dinámica (N · s/m²).

Determine Dp para un tramo horizontal de tubo liso de 0.2 m de longitud, dadas r = 1.23 kg/m³,

m = 1.79 × 10^–5N · s/m², D = 0.005 m, V = 40 m/s, y e = 0.0015 mm. Utilice un método numérico para determinar el factor de fricción.

Obsérvese que los tubos lisos tienen Re < 10⁵, un valor inicial apropiado se obtiene con el uso de la fórmula de Blasius, f = 0.316/Re^0.25. 

(Ejercicio 8.12, tomado del libro Métodos Numéricos para Ingenieros de Steven C. Chapra y Raymond P. Canale. 5Ed. McGraw Hill. Mexico. 2007. Página 218)

SOLUCIÓN

Análisis de la situación

Al leer el problema, vemos que se nos pide calcular la cantidad Dp. Hasta el momento solo tenemos información de L= 0.2 m,  r = 1.23 kg/m³, V = 40 m/s y D=0.0005 m. ¿y cuál es el valor de f, que es el factor de fricción?, muy buena pregunta. Si observamos f está inmersa en la ecuación de Colebrook

 
$$\frac { 1 }{ \sqrt { f }  } =-2{ log }_{ 10 }\left( \frac { \varepsilon  }{ 3.7D } +\frac { 2.51 }{ Re\sqrt { f }  }  \right) $$

En otras palabras, para poder calcular a Dp, primero debemos calcular el factor de fricción. Veamos como determinar a f.

En primer lugar, reemplacemos los datos que tenemos en la ecuación de Colebrook, es decir,

$$\frac { 1 }{ \sqrt { f }  } =-2{ log }_{ 10 }\left( \frac { 0.0015 }{ 3.7*0.005 } +\frac { 2.51 }{ \left( \frac { 1.23*40*0.005 }{ 1.79*{ 10 }^{ -5 } }  \right) \sqrt { f }  }  \right) $$

Al simplificar la ecuación un poco nos quedaría

$$\frac { 1 }{ \sqrt { f }  } =-2{ log }_{ 10 }\left( \frac { 3 }{ 3.7 } +\frac { 4.4929 }{ 24600\sqrt { f }  }  \right) $$

Ahora, para determinar el factor de fricción (f) usaremos el método de Falsa Posición, este método requiere por lo menos una función en una variable y un intervalo que contenga la raíz. La función la construiremos a partir de la ecuación simplificada, es decir,

$$\frac { 1 }{ \sqrt { f }  } =-2{ log }_{ 10 }\left( \frac { 3 }{ 3.7 } +\frac { 4.4929 }{ 24600\sqrt { f }  }  \right) $$

$$\frac { 1 }{ \sqrt { f }  } + 2{ log }_{ 10 }\left( \frac { 3 }{ 3.7 } +\frac { 4.4929 }{ 24600\sqrt { f }  }  \right) = 0$$

$$F(f) = \frac { 1 }{ \sqrt { f }  } + 2{ log }_{ 10 }\left( \frac { 3 }{ 3.7 } +\frac { 4.4929 }{ 24600\sqrt { f }  }  \right)$$

Ahora bien, es claro que el factor de fricción (f) esperado debe satisfacer la igualdad

$$\frac { 1 }{ \sqrt { f }  } =-2{ log }_{ 10 }\left( \frac { 3 }{ 3.7 } +\frac { 4.4929 }{ 24600\sqrt { f }  }  \right) $$

Esto implica que dicho factor al ser evaluado en F(f) esta debe resultar cero, es decir, F(f)=0. Si lo ponemos en otras palabras, el factor de ficción esperado seria aquel valor que sea un intercepto de la gráfica de F(f) con el eje horizontal. Con la gráfica de F(f) y ésta idea en mente, podemos determinar el intervalo que contenga la raíz de F(f) para poder aplicar el método de Falsa Posición. Por comodidad, de ahora en adelante cambiare a ‘f’ por ‘x’, es decir, usare

$$F(x) = \frac { 1 }{ \sqrt { x }  } + 2{ log }_{ 10 }\left( \frac { 3 }{ 3.7 } +\frac { 4.4929 }{ 24600\sqrt { x }  }  \right)$$

Ahora si veamos la gráfica de F(x):

Gráfica 1

En la Gráfica 1, observamos que el intercepto o raíz de F(f) se encuentra en el intervalo [0,1]. Ahora si podemos elegir el intervalo para el método de Falsa Posición, por ejemplo tomemos el intervalo [0.1,1]. Es claro que no podemos incluir el 0, puesto que en este valor la función F(x) tiene una asíntota.

Ya tenemos la función y el intervalo, ahora si apliquemos el método.

Aplicación del método de Falsa Posición

El código en Matlab que implementa el método de Falsa Posición es

f=inline(input('Ingrese funcion:  ','s'));

xa=input('Ingrese extremo inferior del intervalo cerrado: ');

xb=input('Ingrese extremo superior del intervalo cerrado: ');

e=input('Ingrese porcentaje de error: ');

if f(xa)*f(xb)<0

    xrant=xa;   %xrant: raíz antigua

    xr=xa-(f(xa)*(xb-xa)/(f(xb)-f(xa)));   %xr: raíz nueva

    ea=abs((xr-xrant)/xr)*100;

    salir=0;

    while(ea>e && salir==0)       

        if f(xa)*f(xr)<0

            xb=xr;   

        elseif f(xr)*f(xb)<0

            xa=xr;           

        else           

            salir=1;

            xa=xr;

            xb=xr;

        end

        xrant=xr;

        xr=xa-(f(xa)*(xb-xa)/(f(xb)-f(xa)));

        ea=abs((xr-xrant)/xr)*100;

    end

    fprintf('La raiz aproximada es:  %12.10f\n',xr);

else

    fprintf('En el intervalo "No hay Raiz".\n');

end

A este código lo he llamado falsaPosicion.m.

Bueno, ahora si ejecutemos el código e ingresemos la información requerida:

>> falsaPosicion

Ingrese funcion:  1/sqrt(x)+2*log10(3/37+4.4929/(24600*sqrt(x)))

Ingrese extremo inferior del intervalo cerrado: 0.1

Ingrese extremo superior del intervalo cerrado: 1

Ingrese porcentaje de error: 10^(-10)

La raiz aproximada es:  0.2108240043

Por lo tanto, el factor de fricción (f) es 0.2108240043.

Completando la solución

Continuando con la solución de nuestro problema, ya podemos determinar a Dp. Recordemos que de entrada ya teníamos L= 0.2 m,  r = 1.23 kg/m³, V = 40 m/s y D = 0.0005 m. Con el factor de fricción f = 0.2108240043, ya podemos sustituir en

$$\Delta p=f\frac { L\rho { V }^{ 2 } }{ 2D }$$

De modo que

$$\Delta p=0.2108240043*\frac { 0.2*1.23*{ 10 }^{ 2 } }{ 2*0.0005 } $$

$$\Delta p\quad \approx \quad 5186.27\quad Pa$$

Listo, problema solucionado.

Para los amantes de Python, les dejo a continuación el código del método de Falsa Posición:

from sympy import *

funTexto=input('Ingrese funcion: ')

xa=float(input('Ingrese extremo inferior del intervalo cerrado: '))

xb=float(input('Ingrese extremo superior del intervalo cerrado: '))

e=float(input('Ingrese porcentaje de error: '))

f=sympify(funTexto)

if f.subs('x',xa)*f.subs('x',xb)<0:

    xrant=xa   #xrant: raiz antigua

    xr=xa-(f.subs('x',xa)*(xb-xa)/(f.subs('x',xb)-f.subs('x',xa)))   #xr: raiz nueva

    ea=abs((xr-xrant)/xr)*100

    salir=0

    while(ea>e and salir==0):

        if f.subs('x',xa)*f.subs('x',xr)<0:

            xb=xr   

        elif f.subs('x',xr)*f.subs('x',xb)<0:

            xa=xr           

        else:           

            salir=1

            xa=xr

            xb=xr

        xrant=xr

        xr=xa-(f.subs('x',xa)*(xb-xa)/(f.subs('x',xb)-f.subs('x',xa)))

        ea=abs((xr-xrant)/xr)*100

    print('La raiz aproximada es: ',xr)

else:

    print('En el intervalo \"No hay Raiz\".')

Hemos llegado al final de esta entrada, espero les sirva de mucha ayuda.