Matlab - Resolver Iteración para Placa Orificio

Muy buenas a todos,

necesito calcular la relación de diámetros de una placa orificio para medir el caudal de un gas (aire) según la normativa DIN EN ISO 5167.
Para ello se utiliza una ecuación que ha de ser resuelta por iteración. La ecuación es:

C*e*b^2/sqrt(1-b^4) = 4q/(pi*D^2*sqrt(2*p*r);

o escrita de otra manera:

C*e*b^2/sqrt(1-b^4)-4q/(pi*D^2*sqrt(2*p*r)=0;

Donde:

Son incógnitas:
C: Coeficiente de descarga
e: Coeficiente de expansión
b: relación de diámetros d/D

Son datos:
q: caudal
D: diámetro de la tuberia
dp: diferencia de presión entre ambos lados de la placa orificio
rho: densidad del fluido

Para las ecuaciones de coeficiente de descarga y coeficiente de expansión son datos:
p21=p2/p1: relación entre presiones a ambos lados de la placa orificio
k: coeficiente isoentrópico

En realidad solo existe una incógnita "b", ya que "C" y" e" se pueden calcular en función de "b".

Como siempre que he de realizar un cálculo algo más laborioso, he tratado de sacarle partido a las nociones básicas que poseo en MatLab. Desafortunadamente el código que he escrito, a pesar de no dar error, no me sirve. El ordenador comienza a calcular pero tras más de 10 minutos no genera ninguna solución, sigue calculando (se escucha que el ordenador continúa trabajando)..

El código es el siguiente:

>> clear
%Entrada de datos
Re=329546;
p21=137000/139000;
k=1.1378;
q=1.2;
D=0.15;
dp=2000;
rho=0.8031; %fin entrada de datos
syms b; %definicion de variable
e=1-(0.351+0.256*b^4+0.93*b^8)*(1-(p21)^(1/k)); %coef de expansion
C=0.5961+0.0261*b^2-0.216*b^8+0.000521*(10^6*b/Re)^0.7+(0.0188+0.0063*(19000*b/Re)^0.8)*b^(3.5)*(10^6/Re)^(0.3)+(0.043+0.080*exp(-10)-0.123*exp(-7))*(1-0.11*(19000*b/Re)^0.8)*(b^4/(1-b^4))-0.031*((0.94/(1-b)-0.8*(0.94/(1-b))^1.1))*b^(1.3);%coef de descarga
eq=C*e*b^2/sqrt(1-b^4)-4*q/(pi*D^2*sqrt(2*dp*rho));
b=solve(eq,b)

Tal vez esta no sea la mejor manera de solucionar este tipo de problema y exista una función más apropiada.
Estaría muy agradecido si alguien es capaz de ayudarme ya sea con el código que he escrito yo o proponiendo uno totalmente diferente que no use la función solve.

Muchas gracias.

Saludos,

C. Jaimez

Buenas,

Dale tiempo, deja el ordenador calculando, puede tardar mucho. Yo hago iteraciones que me tardan 8h de reloj y solo me hace una iteración de 94 variables, y necesito muuuuchas iteraciones, no te exagero. Que procesador tienes?
Una solucion es monarte un cluster de ordenadores en casa, o pedirlo, o enviar el calculo algun campus y te lo calculen y te lo devuelvan. Ahora estoy probando de ejecutar mi codigo en un cluster en la universidad.

Mi recomendación, dejalo funcionar y ves al trabajo, cuando vuelvas mira si ha acabado.

PD: no mire elcodigo, quizas sea problema de eso, solo te digo que si tarda, es normal, dale tiempo.

David.

Hola David,

gracias por tu comentario. Mi procesador es un i7 y tengo 8GB de Ram. Pero tengo acceso a un cluster en la universidad. No se me pasó por la cabeza que pudiese tardar tanto para este tipo de cálculo, pensé que se trataba de un cálculo más o menos sencillo que el ordenador resolvería deprisa. Algo que no incluí en el código y que puede que acorte bastante el tiempo de cálculo es que la variable a determinar "b" está acotada: 0,10<=b<=0,75
Cómo podría incluir esta restricción en mi código?

Gracias de nuevo.

Christopher

P.D. Si no es mucha molestia igual podrías revisar el código para estar seguro de que es correcto. Gracias!

Podrias usar un método numérico para resolverlo.

Alguna propuesta?

usa método de bisección.

Muchas gracias!