3. Modelado de un Reactor CSTR de Van de Vusse

3. Modelado de un Reactor CSTR de Van de Vusse
5 (100%) 1 vote

Consideremos el siguiente modelo fenomenológico de un reactor de Van de Vusse, al cual se le desea obtener un modelo lineal en diferentes puntos de operación sobre el cual va a trabajar el proceso.

 

Descripción del proceso:

 

El CSTR (reactor perfectamente agitado, del ingles continuous stirred-tank reactor) de van de Vusse, mostrado en la figura a continuación, representa el sistema reacional utilizado para la síntesis de ciclopentenol (representado por B) a partir del ciclopentadieno (representado por A). Este reactor envuelve entonces reacciones en serie y en paralelo entre la entrada del sistema (Reacción A) y el producto obtenido (Reacción B). Sin embargo dado a la fuerte reatividade de los reagentes y produtos, se obtienen dos productos adicionales NO deseados, el diciclopentadieno (D) es produzido como subproduto por la reacción de Diels-Alder y ciclopentenodiol (C) es produzido como um produto consecutivo por la adición de otra molécula de água.

Reactor de Van de Vusse

A\overset{k_1}{\rightarrow}B\overset{k_2}{\rightarrow}C

2A\overset{k_3}{\rightarrow}D

donde los parámetros de la tasa de reacción, vienen dados por:

k_1=5/6 l/min

k_2=5/3 l/min

k_3=1/6 l/min

 

 

Balance General de Masa

Si se asume la densidad y el volumen constante en el interior del reactor tenemos:

\dfrac{dV}{dt}=0 y F\approx F_i

Balance del componente A:

Interior del reactor = Flujo de entrada – Flujo de Salida – sale reaccion 1 – sale reacción 2

\dfrac{d(VC_A)}{dt}=F(C_{Af}-C_A)-Vk_1C_A-VK_3C_A^2

Y debido a que consideramos el volumen constante, el balance de masa del componente A se resume a:

\dfrac{d(C_A)}{dt}=\dfrac{F}{V}(C_{Af}-C_A)-k_1C_A-K_3C_A^2

Balance del componente B:

De la misma forma podemos obtener el balance de masa para el componente B

\dfrac{d(C_B)}{dt}=-\dfrac{F}{V}C_B+k_1C_A-K_2C_B

Balance del componente C:

\dfrac{d(C_C)}{dt}=-\dfrac{F}{V}C_C+k_2C_B

Balance del componente D:

\dfrac{d(C_D)}{dt}=-\dfrac{F}{V}C_D+\dfrac{1}{2}k_3C_A^2

Como deseamos obtener el modelo del Componente B con relación al producto de entrada Componente A, observemos que estos dos balances de energia NO dependen de los balances C y D, por lo tanto, solo nos van a interesar estas dos ecuaciones diferenciales:

\dfrac{d(C_A)}{dt}=\dfrac{F}{V}(C_{Af}-C_A)-k_1C_A-K_3C_A^2

\dfrac{d(C_B)}{dt}=-\dfrac{F}{V}C_B+k_1C_A-K_2C_B

Estas dos ecuaciones representan el comportamiento fenomenológico no lineal de nuestro reactor de Van de Vusse. Ahora como deseamos linealizar este sistema en torno a un punto de equilibrio o punto de operación el próximo paso sera encontrar dicho punto. Para lograr esto hacemos todas las derivadas (cambios) iguales a cero, lo que significa que el sistema no está cambiando y se encuentra en un estado estacionario.

A partir de la ecuación de C_A obtenemos lo siguiente:

0=\dfrac{F_s}{V}(C_{Afs}-C_{As})-k_1C_{As}-K_3C_{As}^2

-K_3C_{As}^2+\left(-k_1-\dfrac{F_s}{V}\right)C_{As} +\dfrac{F_s}{V}C_{Afs}=0

Resolviendo la ecuación cuadrática anterior utilizando la formula general y escogiendo la raiz positiva, tenemos que:

C_{As}=\dfrac{-\left(k_1+\dfrac{F_s}{V}\right)}{2k_3}+\dfrac{\sqrt{\left(k_1+\dfrac{F_s}{V}\right)^2+4k_3\dfrac{F_s}{V}}}{2k_3}

Y para el componente B tenemos:

0=-\dfrac{F_s}{V}C_{Bs}+k_1C_{As}-K_2C_{Bs}

C_{Bs}=\dfrac{k_1C_{As}}{\dfrac{F_s}{V}+K_2}

La entrada de nuestro sistema es la tasa de disolución F/V l/min. Si hacemos un barrido de nuestra entrada desde 1 1/min hasta 10 1/min. Obtendremos las curvas que relacionan la entrada y la salida del sistema en su estado estacionario.

Van de Vusse Estacionario
Si observamos el comportamiento del Componente B, vemos que tiene una inversión de fase, que con una entrada de 1.292 1/min se obtiene la máxima producción del componente B (Lo más deseado) y que con una entrada superior la producción de B comienza a caer.

Linealización del sistema de ecuaciones no lineales en Espacio de Estados utilizando el Jacobiano

Vamos a obtener una representación en espacio de estados de nuestro reactor.

\dot{x}=Ax+Bu

y=Cx

Definimos cuales serán nuestras variables de estado:

x_1=C_A

\dot{x}_1=\dfrac{dC_A}{dt}

x_2=C_B

\dot{x}_2=\dfrac{dC_B}{dt}

u= F/V

Si renombramos las dos ecuaciones de balance de masa como función 1 (f_1) y función 2 (f_2).

\dfrac{d(C_A)}{dt}=f_1(C_A,C_B,F/V)=\dfrac{F}{V}(C_{Af}-C_A)-k_1C_A-K_3C_A^2

\dfrac{d(C_B)}{dt}=f_2(C_A,C_B,F/V)=-\dfrac{F}{V}C_B+k_1C_A-K_2C_B

Jacobiano:

A=\begin{bmatrix}  \dfrac{\partial f_1}{\partial x_1} & \dfrac{\partial f_1}{\partial x_2}\\ \\  \dfrac{\partial f_2}{\partial x_1} & \dfrac{\partial f_2}{\partial x_2}  \end{bmatrix},B=\begin{bmatrix}  \dfrac{\partial f_1}{\partial u} \\ \\  \dfrac{\partial f_2}{\partial u}  \end{bmatrix}

Matriz A:

A_{11}=\dfrac{\partial f_1}{\partial x_1}=-\dfrac{F_s}{V}-k_1-2k_3C_{As}

A_{12}=\dfrac{\partial f_1}{\partial x_2}=0

A_{21}=\dfrac{\partial f_2}{\partial x_1}=k_1

A_{22}=\dfrac{\partial f_2}{\partial x_2}=-\dfrac{F_s}{V}-k_2

Matriz B:

B_{11}=\dfrac{\partial f_1}{\partial u}=C_{Afs}-C_{As}

B_{21}=\dfrac{\partial f_2}{\partial u}=-C_{Bs}

El modelo en espacio de estado viene dado por:

\dot{x}=\begin{bmatrix}  -\dfrac{F_s}{V}-k_1-2k_3C_{As} &0\\ \\  k_1 & -\dfrac{F_s}{V}-k_2  \end{bmatrix}x+\begin{bmatrix}  C_{Afs}-C_{As} \\ \\  -C_{Bs}  \end{bmatrix}u

y=\begin{bmatrix}  0&1  \end{bmatrix}x

Teniendo la representación en el espacio de estado de nuestro sistema, vamos a seleccionar diferentes puntos de operación para analizar el comportamiento de este proceso. Veremos tres casos diferentes. En el caso 1 veremos como aparece una respuesta inversa o respuesta de fase no mínima en la concentración de B mientras que en el caso 2 esta respuesta inversa no aparece. Y en el caso 3 se ilustra el punto de operación donde se alcanza la máxima producción del componente B.

CASO 1

Para este punto de operación vamos a considerar la tasa de flujo \dfrac{F_s}{V}=\dfrac{4}{7}min^{-1}. Notemos que esta tasa de flujo se ubica en el lado izquierdo del máximo pico del gráfico de estado estacionario que relaciona la entrada salida.

Inicialmente encontramos el estado estable del componente A (C_{As}) y el componente B(C_{Bs}).

C_{As}=\dfrac{-\left(k_1+\dfrac{F_s}{V}\right)}{2k_3}+\dfrac{\sqrt{\left(k_1+\dfrac{F_s}{V}\right)^2+4k_3\dfrac{F_s}{V}}}{2k_3}

C_{As}=\dfrac{-\left(5/6+4/7\right)}{2(1/6)}+\dfrac{\sqrt{\left(5/6+4/7\right)^2+4(1/6)(4/7)}}{2(1/6)}

C_{As}=3\dfrac{mol}{litro}

C_{Bs}=\dfrac{k_1C_{As}}{\dfrac{F_s}{V}+K_2}

C_{Bs}=\dfrac{3(5/6)}{4/7+5/3}

C_{Bs}=1.117\dfrac{mol}{litro}

Con esto obtenemos la siguiente representación en el espacio de estado

\dot{x}=\begin{bmatrix}  -2.4048 &0\\ \\  0.8333 & -2.2381  \end{bmatrix}x+\begin{bmatrix}  7 \\ \\  -1.117  \end{bmatrix}u

y=\begin{bmatrix}  0&1  \end{bmatrix}x

Y podemos obtener la función de transferencia de la misma forma que lo hemos visto en el curso de sistemas dinámicos lineales usando G(s)=C(sI-A)^{-1}B.

G(s)=\dfrac{-1.117s+3.1472}{s^2+4.6429s+5.3821}

Si notamos, vemos que esta función de transferencia tiene un cero positivo, lo que provoca la respuesta inversa o respuesta de fase no mínima.

Cero=\dfrac{3.1472}{1.117}=2.8175

Al final del post, puedes encontrar el código de implementación en Matlab para que lo copies y lo pegues en tu computador.

 

Veamos una comparación de los dos sistemas en el punto de operación linealizado

Van de Vusse caso 1

Ahora una comparación alejándonos más del punto de operación linealizado. Veamos como el sistema lineal comienza a dejar de valer.

Van de Vusse caso 1

CASO 2

Para este punto de operación vamos a considerar la tasa de flujo \dfrac{F_s}{V}=2.8744min^{-1}. Notemos que esta tasa de flujo se ubica en el lado derecho del máximo pico del grafico de estado estacionario que relaciona la entrada salida.

Procedemos a encontrar el punto estable de ambos componentes de la misma forma que lo hicimos en el caso 1:

C_{As}=6.0870\dfrac{mol}{litro}

C_{Bs}=1.117\dfrac{mol}{litro}

Con esto obtenemos la siguiente representación en el espacio de estado

\dot{x}=\begin{bmatrix}  -5.7367 &0\\ \\  0.8333 & -4.5411  \end{bmatrix}x+\begin{bmatrix}  3.9130 \\ \\  -1.117  \end{bmatrix}u

y=\begin{bmatrix}  0&1  \end{bmatrix}x

Y podemos obtener la función de transferencia de la misma forma que lo hemos visto en el curso de sistemas dinámicos lineales usando G(s)=C(sI-A)^{-1}B.

G(s)=\dfrac{-1.117s-3.1472}{s^2+10.2778s+26.0508}

En este caso ya se tiene el cero en el semiplano negativo, por lo tanto la respuesta inversa desaparece.

Cero=-\dfrac{3.1472}{1.117}=-2.8175

Van de Vusse caso 2

CASO 3

El punto de operación óptimo para la producción del componente B se logra utilizando la siguiente tasa de disolución:

\dfrac{F_s}{V}=1.2921min^{-1}

C_{As}=4.4949\dfrac{mol}{litro}

C_{Bs}=1.266\dfrac{mol}{litro}

Con esto obtenemos la siguiente representación en el espacio de estado

\dot{x}=\begin{bmatrix}  -3.6337 &0\\ \\  0.8333 & -2.9588  \end{bmatrix}x+\begin{bmatrix}  5.5051 \\ \\  -1.2660  \end{bmatrix}u

y=\begin{bmatrix}  0&1  \end{bmatrix}x

Y podemos obtener la función de transferencia de la misma forma que lo hemos visto en el curso de sistemas dinámicos lineales usando G(s)=C(sI-A)^{-1}B.

G(s)=\dfrac{-1.2660s}{s^2+6.5825s+10.7217}

El cero de esta función de transferencia se ubica en el origen. En este caso la condición de estado estable no cambia en la concentración de B para pequeños cambios en la entrada. Esto es debido a la ganancia cero de la función de transferencia.

 

Van de Vusse caso 3

Bibliografica

  • Bequette, B. (1998). Process Dynamics: Modeling, Analysis, and Simulation. New Jersey: Prentice Hall

Código de Implementación

MATLAB:

SIMULINK:

A continuación se presenta el código de implementación para que lo puedas copiar y pegar en tu computador. Recuerda que para poder ver el código, debes compartir el contenido de esta entrada con alguno de los tres botones mostrados abajo. Asi ayudas a que más personas conozcan la página y se beneficien con esta información.

>>>>> Descargar Archivos de Simulación <<<<<<

Principal:

 

 

Función del Modelo:

 

Si te sirvió de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.




Esta entrada tiene 2 comentarios

  1. oi Sérgio td bem? fiz a disciplina do Maurício com vc esse semestre na Ufrj mas aluna da puc. Eu estava conferindo seu código com o meu pois tive alguns erros nos meus trabalhos e queria entender melhor algo: no seu código, o k1, k2 e k3 são constantes, certo? isso é, não variam com a temperatura como no exercício da disciplina, certo? ou esses k1, k2 e k3 do seu código se referem aos k1, k2 e k3 iniciais?

    1. Olá Sarah, tudo certo!!. É correto teu analise. Neste exemplo aqui, são consideradas as constantes cineticas K1, K2, e K3 como constantes que são aqualas que aparecem no código e não variam com a Temperatura. O modelo do exemplo desta entrada, não está conseiderando-se o balaço de energia da jaqueta.
      Na outra entrada, está o reator de VV para o caso MIMO, que é igual ao visto na disciplina, ali sim consideramos a jaqueta e as constantes variam com relação à temperatura. Podes dar click aqui, caso tu não tenhas visto o outro exemplo. Qualquer dúvida é só me procurar. Grande abraço.

Deja un comentario

Menú de cierre