Mecánica Computacional Vol XXV, pp. 1313-1334 Alberto Cardona, Norberto Nigro, Victorio Sonzogni, Mario Storti. (Eds.)
Santa Fe, Argentina, Noviembre 2006
MODELIZACIÓN NUMÉRICA DE UN COMPRESOR MONOCILINDRICO DE DESPLAZAMIENTO POSITIVO
G. Bourgesa,1, J. Eliacha,2, M. Medinab
a Escuela de Ingeniería Mecánica. Fac. de Ciencias Exactas, Ingeniería y Agrimensura. UNR.
Beruti 2109, (2000) Rosario, Argentina. gbourges@fceia.unr.edu.ar1, jeliach@gmail.com2.
b Consejo de Investigaciones. Universidad Nacional de Rosario. Escuela de Formación Básica. Fac. de Ciencias Exactas, Ingeniería y Agrimensura. UNR. Pellegrini 250. (2000) Rosario. Argentina.
mmedina@fceia.unr.edu.ar. Palabras claves: compresor monocilíndrico, metano, simulación, MATLAB
Resumen. En el presente trabajo se simula el proceso de compresión de gas metano (CH4) por medio de un compresor de desplazamiento positivo a pistón (movimiento alternativo), analizando el comportamiento del gas. Se modelizó el sistema compresor–gas modulándolo por medio de componentes. Se programó el simulador en entorno MATLAB, utilizando las herramientas numéricas de dicho software.
Se muestra mediante un ejemplo numérico la aplicación del simulador, con determinadas condiciones de contorno, para predecir el comportamiento del gas. Los resultados obtenidos muestran que el modelo es capaz de predecir el comportamiento publicado en la bibliografía. Este simulador puede resultar muy útil durante el diseño de las partes del compresor porque permite la experimentación numérica con diferentes configuraciones para evaluar la incidencia de cada componente y conseguir la combinación óptima de los mismos.
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
1314
G. BOURGES, J. ELIACH, M. MEDINA
1 INTRODUCCION
Los compresores son máquinas que tienen por finalidad aportar una energía a los fluidos compresibles (gases y vapores) sobre los que operan, para hacerlos fluir aumentando al mismo tiempo su presión. En esta última característica precisamente, se distinguen de las soplantes y ventiladores que manejan grandes cantidades de fluidos compresibles (aire por ejemplo) sin modificar sensiblemente su presión, con funciones similares a las bombas de fluidos incompresibles (Fernández Diéz, 2002) .
El gas metano (principal componente del Gas Natural Comprimido ó G.N.C.) es uno de los combustibles mas requeridos dentro de la República Argentina. El transporte de este gas se realiza mayormente por medio de grandes gasoductos, que van desde los lugares de extracción del gas hacia los distintos centros de consumo. Dichos gasoductos conducen enormes volúmenes de gas y requieren estaciones compresoras y de recompresión para mantener el flujo constante necesario en las condiciones de presión adecuadas (Faires, 1999) .
Los compresores son de dos tipos generales: de movimiento alternativo (o de cilindro y pistón), y de movimiento rotativo. Estos últimos pueden ser de desplazamiento positivo en cuyo caso el volumen entregado es constante (por ejemplo compresores a tornillo, a paletas, etc.) o centrífugo en donde el volumen es variable en función de la presión. Los de desplazamiento positivo se utilizan generalmente cuando se requieren elevadas relaciones de presión (presión de salida del compresor respecto de la presión de ingreso) (Faires, 1999). En este trabajo se analiza un compresor de desplazamiento positivo del tipo cilindro y pistón.
1.1 Descripción y funcionamiento
1.1.1 Descripción del equipo
El compresor simulado está compuesto por los siguientes elementos (Figura 1): el cilindro de compresión; el pistón ó embolo de compresión, acoplado al sistema biela-cigüeñal; y la tapa del cilindro en donde se ubican las válvulas y conductos de admisión y escape del gas. El cilindro posee en su superficie exterior aletas de enfriamiento las cuales tiene como finalidad aumentar la superficie de intercambio de calor entre el cilindro y el medio ambiente.
El movimiento angular del cigüeñal es transformado en movimiento alternativo por el sistema cigüeñal-biela-pistón. En su recorrido, el pistón posee dos posiciones extremas: Punto Muerto Superior ó P.M.S. (Figura 1) que corresponde al menor volumen, y el Punto Muerto Inferior ó P.M.I. (Figura 1) correspondiente a un mayor volumen dentro del cilindro. El menor volumen encerrado por el cilindro se lo nombra como volumen ó espacio nocivo, estimándose entre un 3% a un 10% del volumen desplazado según el modelo del compresor (Fernández Diéz, 2002). El volumen desplazado por el pistón durante su recorrido corresponde a la diferencia entre el volumen total ó máximo del cilindro y el volumen nocivo.
La existencia de un volumen nocivo provoca un retraso en el ingreso de gas durante el proceso de aspiración debido a que el gas almacenado en dicho volumen a la presión Pescape debe expandirse hasta la presión Padmisión (Figura 3) antes de permitir la entrada de gas al cilindro. Sin embargo, su efecto es doble dado que si bien por un lado disminuye el volumen de aspiración, por otro lado ahorra energía, ya que la expansión produce un efecto motor sobre el pistón; es por ello que se puede considerar que ambos efectos se compensan desde el punto de vista energético (Fernández Diéz, 2002).Las válvulas de admisión y escape (Figura 1) permitirán el ingreso ó egreso del gas respectivamente, según como sea la relación de presiones entre la cámara de compresión y el interior de los respectivos conductos de admisión ó escape.
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
Mecánica Computacional Vol XXV, pp. 1313-1334 (2006)
1315
Conducto de Admisión
Válvula de
ALETAS DE Admisión ENFRIAMIENTO
P.M.S
Desplazamiento pistón
P.M.I.
BIELA
Conducto de Escape
Válvula de Escape
Volumen Nocivo
CILINDRO
PISTÓN
Giro del cigüeñal
CIGÜEÑAL
Figura 1: Esquema del compresor.
1.1.2. Funcionamiento y ciclo de compresión
En un compresor de gas de movimiento alternativo el fluido ingresa dentro de la cámara de compresión a una presión Padmisión, descargándolo al final del proceso a una presión Pescape superior (Figura 3) (Fernández Diéz, 2002). El proceso de compresión en un compresor a émbolo consta de cuatro etapas, las cuales se describen a continuación:
1º Ingreso de gas al cilindro: ocurre desde el momento en el cual se produce la apertura de la válvula de admisión ó A.V.A (punto 1). Este proceso se prolonga hasta su cierre (punto 2) y se denota con las siglas C.V.A. (Cierre de Válvula de Admisión), siendo la presión dentro del cilindro durante este período constante e igual a la presión en el conducto de admisión. En proceso reales, la A.V.A. se retrasará hasta que la presión dentro del cilindro sea (PadmisiónΔPadmisión), es decir una presión menor a la existente en el conducto de admisión. La diferencia ΔPadmisión corresponde al hecho de que el gas deberá vencer tanto la tensión del resorte de la válvula de admisión como la resistencia al paso de flujo debido al estrangulamiento que ella provoca. Por lo tanto el gas sufrirá una expansión ΔPadmisión en su paso por el orificio de la válvula de admisión. Se supone que la fuente de ingreso de gas es lo suficientemente grande para no producir efectos pulsantes.
El C.V.A. coincide, en un proceso teórico, con una de las posiciones extremas del pistón, el Punto Muerto Inferior ó P.M.I. En un proceso real, el cierre de la válvula de admisión (C.V.A.) se produce instantes previos a la llegada del pistón al P.M.I.
2º Proceso de compresión: Esta etapa inicia cuando el pistón se encuentra en el P.M.I. (punto 2) desplazándose hacia el P.M.S., y se prolonga hasta la apertura de la válvula de escape ó A.V.E. (punto 3).
El proceso de compresión del gas es descrito por la siguiente ecuación:
pV k = cte
(1)
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
1316
G. BOURGES, J. ELIACH, M. MEDINA
Cuando el valor de k es la unidad ( k = 1) la compresión es isotérmica, es decir el gas se
comprime a temperatura constante, siendo la curva de compresión como la curva 1-2 (Figura
2) . La ecuación (1) queda reducida a:
pV = cte
(2)
Cuando el valor de k coincide con la constante γ = Cp/Cv (relación entre los calores específicos del gas), la compresión es adiabática siendo su curva la 1-2’ (Figura 2). Es evidente que, en este último caso, la cantidad de trabajo necesaria para comprimir y descargar el gas es mayor que cuando la compresión es isotérmica. En consecuencia, el trabajo necesario para comprimir un gas disminuye a medida que se reduce el valor del exponente k (Faires, 1999).
En un caso real, la curva de compresión se encuentra entre las dos curvas límites (1-2 y 12’), encontrándose el valor del coeficiente k dentro los límites 1 < k < γ . En estos casos se dice
que la compresión es politrópica. Nótese que durante una compresión isotérmica se extrae todo el calor equivalente al
trabajo de compresión, y el gas fluye entonces del compresor con la misma energía interna que poseía al ingresar al mismo. En la compresión adiabática no se transmite calor al exterior, y el gas fluye con un aumento de energía interna equivalente al de compresión. Finalmente, en la compresión politrópica existe una cierta absorción de calor por parte del medio exterior, y un determinado incremento en la energía interna y en la temperatura (Faires, 1999). La disipación del calor proveniente del aumento de energía interna que produce el trabajo de compresión puede obtenerse por medio de sistemas de circulación de agua ó aire de enfriamiento alrededor del cilindro de compresión.
Es importante señalar que no hay ninguna ventaja en que el gas comprimido salga del compresor con una energía interna mayor a la que se obtendría con una compresión isotérmica (Faires, 1999).
Presión
2
2'
pV = c
pVγ= C
1
Volumen
Figura 2: Comparación del trabajo entre una compresión isotérmica y una compresión adiabática.
3º Egreso del volumen de gas comprimido. Esta etapa se produce mientras la válvula de escape permanece abierta, es decir entre A.V.E. (punto 3) y C.V.E. (punto 4) En un caso ideal, el C.V.E. coincidirá con el final de carrera del émbolo de compresión ó P.M.S. (Punto Muerto Superior). También, en un caso ideal, la presión dentro del cilindro en esta etapa será constante (Pescape), manteniéndose hasta el cierre de la válvula de escape (Figura 3). En un caso real la presión variará durante el proceso de escape del gas, siendo siempre superior a la presión del gas en la tubería de escape (Figura 4). Esa diferencia de presión ΔPescape corresponde a los mismos factores que disminuye la presión dentro del cilindro durante la
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
Mecánica Computacional Vol XXV, pp. 1313-1334 (2006)
1317
admisión, solo que ahora la presión dentro del cilindro debe ser superior al la de la tubería de escape para que exista flujo de gas.
4º Expansión sin intercambio de masa. Este proceso comienza una vez cerrada la válvula de escape ó C.V.E. (punto 4) y se prolonga hasta la apertura de la válvula de admisión ó A.V.A. (punto 1).
Pescape
Proceso de escape
4
3
P.M.I.=C.V.A.
P.M.S. = C.V.E A.V.A. A.V.E.
Padmisión
Proceso de admisión 1
2
Desplazamiento de pistón
Figura 3: Ciclo de trabajo teórico de un compresor.
Presión interna
Pescape
Proceso de escape
Padmisión
Proceso de admisión
Desplazamiento de pistón
Figura 4: Comparación entre ciclo de trabajo teórico (curvas y rayado verde) y un ciclo de trabajo real (curvas y rayado rojo) en un compresor.
1.2. Simulador
El objetivo de este trabajo fue desarrollar un software que permita evaluar las variables presión y temperatura de gas metano (CH4), y la variación de estos parámetros durante los procesos de admisión, compresión y escape del mismo en un compresor de desplazamiento positivo a pistón. Para ello se incluyó en el código la mayor cantidad de parámetros considerados relevante por los autores y la bibliografía (Faires, 1999; Jones, 1997; Fernández Diéz, 2002).
El software desarrollado incorpora al cálculo de presión y temperatura del gas encerrado en el cilindro, una ley de variación volumétrica del espacio ocupado por el gas, un modelo de dinámica de apertura y cierre de las válvulas de admisión y escape de gas, y un modelo de pérdidas térmicas entre el gas y el cilindro de compresión y entre este último y el medio exterior.
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
1318
G. BOURGES, J. ELIACH, M. MEDINA
1.2.1. Datos requeridos por el simulador
Los datos requeridos por el simulador son: Características del gas: Condiciones del gas en la admisión y escape (presión y temperatura). Constantes típicas del gas: calor específico a presión constante (Cp) y a volumen constante (Cv), constante del gas (R). Geometría del sistema: Diámetro del cilindro, volumen nocivo del cilindro, espesor de la pared del cilindro, longitud de biela, longitud de cigüeñal, diámetro de válvulas, espesor y altura de los aros de pistón, cantidad de aros de pistón y superficie total de aletas de enfriamiento. Otros datos: Características del material constructivo del cilindro (conductividad térmica, calor específico y densidad). Características de las válvulas y tapa de cilindro y condiciones iniciales del gas encerrado en el cilindro (temperatura inicial y masa inicial encerrada).
1.2.2. Datos obtenidos del simulador
Los datos obtenidos del simulador son: la presión, la temperatura y densidad que posee el gas encerrado dentro del cilindro durante el período de tiempo simulado. Otro de los datos obtenidos es el volumen encerrado por el cilindro en cada intervalo de tiempo, dato utilizado para calcular la densidad del gas.
El desarrollo del simulador constó de dos etapas. En la primera se realizaron modelos matemáticos de los componentes que forman el compresor simulado. En la segunda etapa se programó el simulador en el entorno MATLAB®, utilizando las herramientas numéricas de dicho software.
2 MODELADO
En esta sección se presenta el modelado del compresor. En la Figura 5 se muestra un esquema del compresor donde se observan las zonas a tener en cuenta para la modelización.
Propiedades del gas
V.A.
Modelo de movimiento V.E. de válvulas
P.M.S
Modelo de pérdidas térmicas
P.M.I.
Modelo de movimiento del pistón
Figura 5: Modelos desarrollados.
Teniendo como objetivo el cálculo de las propiedades del gas (densidad, presión y temperatura) se plantea un balance energético a través de un volumen de control ubicado en el espacio ocupado por el gas (Figura 6), que coincide con el volumen desplazado por el pistón mas el volumen nocivo, se tiene la siguiente ecuación:
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
Mecánica Computacional Vol XXV, pp. 1313-1334 (2006)
1319
•
••
•
Eint + W + Q pared + Q aros + ΔH = 0
(3)
Donde cada término de la ec. (3) representa:
•
Eint
=
d (Eint ) dt
Variación de energía interna del gas respecto del tiempo.
•
W
=
•
Pcil (t) ⋅Vcil
=
Pcil
(t
)
⋅
d
(Vcil dt
(t
))
Potencia transferida por el pistón al gas en forma de
trabajo. Se adopta signo positivo cuando el volumen disminuye (entrega energía al gas).
•
Q pared Potencia (cantidad de calor por unidad de tiempo) disipada entre el gas y el medio
exterior, a través de la pared del cilindro.
•
Qaros Potencia disipada en forma de calor por unidad de tiempo debido al rozamiento
entre los aros de pistón y la superficie interna del cilindro.
•
ΔH = m h Variación de la entalpía, entrante (válvula de admisión) ó saliente (válvula de
escape) a través de las válvulas.
h
Pcil(t), Vcil(t), Tcil(t)
Volumen de control
Qpared
Qaros
W Qaros
Qpared
Desplazamiento pistón
Figura 6: Balance energético.
La ecuación (3) será discretizada en el de tiempo, resolviéndose numéricamente en cada paso de tiempo los valores de energía interna Eint y de masa del gas encerrado en el cilindro. A partir de estos valores se calculan las propiedades del gas (presión, densidad y temperatura) del siguiente modo:
Cálculo de la densidad: Cálculo de temperatura: Cálculo de presión:
ρ gas
=
mgas Vcil
⎡ Kg ⎤ ⎢⎣ m3 ⎦⎥
[ ] ΔTcil
=
ΔEint Cv
ºK
[ ] Pcil = ρgas.R.Tcil Pa
(4) (5) (6)
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
1320
G. BOURGES, J. ELIACH, M. MEDINA
En las secciones siguientes se describen los modelos utilizados para el cálculo de los términos de la ecuación (3).
•
2.1 Cálculo de W . Potencia transferida por el pistón.
•
En este ítem se plantea el modelo matemático de W a través del cálculo de la variación
temporal
del
volumen
•
Vcil
=
d (Vcil (t)) dt
.
Para
el
valor
de
Pcil (t) se
adopta
el
obtenido
en
el
paso
de tiempo anterior.
•
Para el cálculo de Vcil es necesario calcular el desplazamiento del pistón de compresión.
Dicho movimiento se realiza por medio de un sistema tipo biela manivela. Dado que el
cálculo de las propiedades del gas, así como el modelo de pérdidas térmicas requieren
información acerca de la geometría del cilindro, es de vital importancia conocer la posición
del pistón en cada instante. Se sabe además que la posición del pistón esta directamente
relacionada con el ángulo girado por el cigüeñal (Heywood, 1988; Nigro, 1999).
Considerando un esquema del compresor como el mostrado en la Figura 7, donde D es el diámetro del cilindro, L la longitud de la biela, a el radio del cigüeñal,φ = φ(t) es el ángulo
girado por el cigüeñal en cada instante y P.M.I. (punto muerto inferior) y P.M.S. (punto
muerto superior) son las posiciones extremas que toma el pistón en cada ciclo.
D
P.M.S
l Δx Dtp
P.M.I.
L
φ a
Figura 7: Esquema del compresor.
De la bibliografía (Heywood, 1988; Nigro, 1999) se sabe que la distancia entre el eje del cigüeñal y el perno del pistón en cada instante será:
s(t) = a.cos(φ(t)) + L2 − a2.sen(φ(t))2
(7)
De (7) se ve que el menor valor de s se dará cuando φ = π , es decir cuando el pistón se
encuentre en el P.M.I., mientras que el máximo será cuando φ = 0 , que ocurre cuando el
pistón se encuentre en el P.M.S. En otras palabras, la distancia s varía entre los valores extremos: (L − a) ≤ s ≤ (L + a) , por lo que la distancia a ser recorrida por el pistón entre el
P.M.I. y el P.M.S. es de
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
Mecánica Computacional Vol XXV, pp. 1313-1334 (2006)
1321
Δx = (L + a) − (L − a) = 2.a
(8)
Si lo que se desea conocer es la distancia entre la tapa de válvulas y la cabeza del pistón, entonces debe realizarse la siguiente operación:
Dtp (t) = l + (L + 2.a) − s(t)
(9)
Esta distancia varía entre l ≤ Dtp ≤ l + 2.a .
La tasa de variación de la distancia s respecto al ángulo girado por el cigüeñal
ds dφ
es
ds dφ
=
−⎛⎝⎜⎜ a.sen(φ(t))
+
2.
a2.sen(2.φ(t)) ⎟⎞ L2 − a2.sen(φ(t))2 ⎠⎟
(10)
De (7), (9), (10) se observa que la única variable independiente es φ = φ(t) , es decir el
ángulo girado por el cigüeñal. Conociendo la velocidad de rotación del motor ω, y sabiendo que ω = dφ (siendo t la variable tiempo) es constante, puede entonces obtenerse el valor de
dt φ para cada instante de tiempo, teniendo en cuenta las condiciones iniciales:
dφ = ω.dt
(11)
De (10) y (11) podemos obtener la tasa de variación de la distancia s respecto al tiempo t, es decir la velocidad instantánea del pistón:
•
s
=
ds dt
=
−ω.⎜⎛⎜⎝ a.sen(φ(t))
+
2.
a2.sen(2.φ(t)) ⎞⎟ L2 − a2.sen(φ(t))2 ⎟⎠
(12)
El volumen encerrado en el cilindro en cada instante es:
Vcil
(t)
=
Vc
+
π .D2 4
(L
+
a
−
s(t))
(13)
Por lo tanto la tasa de variación de volumen en respecto de φ es:
d (V cil (t)) dφ (t )
=
⎛⎜⎜⎝
π
.D 4
2
⎞⎟⎠⎟
ds dφ
=
⎛⎜⎝⎜
π
.D 4
2
⎟⎞⎠⎟
ds dt
1 ω
(14)
Y la tasa de variación de volumen en respecto del tiempo t es:
•
V
=
d (V cil (t)) dt
=
⎛⎜⎝⎜
π
.D 4
2
⎟⎟⎠⎞
ds dt
(15)
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
1322
G. BOURGES, J. ELIACH, M. MEDINA
•
2.2 Cálculo de Q pared . Modelo de pérdidas térmicas
De lo visto en la ecuación de balance energético (3), se desprende la necesidad de conocer el flujo de calor por unidad de tiempo que se disipa a través de la pared del cilindro. Para ello se plantea un modelo, en el cual se tiene en cuenta la transferencia de calor entre la masa de gas encerrada dentro del cilindro de compresión y la pared del cilindro, y entre esta última y el medio exterior. No se considera aquí el intercambio de calor a través de la tapa del cilindro ni el que pudiera existir a través del pistón. Esta simplificación se basa en el hecho de suponer que, como el compresor simulado posee aletas de enfriamiento solo en su superficie cilíndrica exterior, es a través de dicha superficie que se dará el mayor flujo de calor entre el medio exterior y el gas, y por lo tanto podemos despreciar la contribución de las otras superficies. Dichas aletas de enfriamiento son tenidas en cuenta en el cálculo solo por su superficie de intercambio de calor con el aire. También incluye el aporte de energía proveniente de la fricción entre los aros de pistón y la pared interna del cilindro.
La ecuación que describe el flujo de calor entre el la pared interior del cilindro y el gas encerrado en el mismo, es la siguiente:
( ) •
Q pared = Ωint\_ cil .hc . Tcil − Tpared \_ int
(16)
Donde Ωint\_cil es la superficie interior del cilindro en contacto con el gas, hc el coeficiente
pelicular de transferencia de calor por convección entre el gas y la superficie, Tcil y Tpared \_int la temperatura del gas y de la superficie interna. La superficie Ωint\_cil depende de la posición del pistón en cada instante, pues es zona de contacto entre el gas y la pared interna del cilindro, y tiene la forma:
Ωint\_ cil = π .D.Dtp (t)
(17)
Donde D es el diámetro del cilindro y Dtp(t) es la distancia entre la tapa de válvulas y la cabeza del pistón vista en la ec. (9).
De (17) se observa que conociendo la superficie de intercambio Ωint\_cil , es necesario
( ) encontrar el valor del coeficiente pelicular hc y del gradiente térmico entre el gas y la pared
T − T cil
pared \_ int
.
Respecto
a
esta
expresión,
se
debe
calcular
también
la
temperatura
de
la
superficie interna del cilindro Tpared \_int , pues la temperatura del gas Tcil es conocida de (5) del
paso de tiempo anterior.
Con el objetivo de obtener mayor claridad se desarrollan en los dos ítem. siguientes el
cálculo del coeficiente pelicular hc y de la temperatura de pared Tpared \_int .
2.2.1 Cálculo de hc entre el gas y la superficie interior del cilindro
En este ítem, se analiza el cálculo del coeficiente pelicular hc de transferencia de calor entre el gas y la superficie interior del cilindro. Se supone que el intercambio de calor entre el gas y la pared interior del cilindro se realiza principalmente por convección, despreciándose los efectos de la conducción de calor y de radiación entre el gas y el cilindro. También se
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
Mecánica Computacional Vol XXV, pp. 1313-1334 (2006)
1323
supone que el gas en estudio permanece siempre como una mezcla uniforme con iguales propiedades (densidad, presión y temperatura) cualquier punto del volumen de control.
El coeficiente pelicular de transmisión de calor entre el gas y la pared hc se obtiene entonces de la siguiente forma ( Heywood, 1988; Faires, 1999; Kreith, 2001; Bird, 2001):
hc\_gas
=
0.035. Re 0.8 D
.Pr0.33 .k
(18)
Siendo:
Re = S pistón ⋅ D
(19)
ν
La ecuación (19) representa el numero adimensional de Reynolds, donde S pistón es la
velocidad promedio del pistón (Heywood, 1988; Nigro, 1999). Para el número adimensional de Prandtl se adopto el valor Pr = 0,72 ( Heywood, 1988; Faires, 1999; Kreith, 2001; Bird,
2001).
K gas
=
ν.C p Pr
(20)
2.2.2 Calculo de la temperatura de la superficie interna del cilindro Tpared \_int
Este ítem tiene como objetivo la obtención de la temperatura de la superficie interna del cilindro T , pared \_int necesaria para conocer el flujo de calor entre el gas y el exterior. Dicho
flujo dependerá tanto de las condiciones internas del cilindro (temperatura del gas (Tcil ), área de contacto entre el gas y el cilindro ( Ωint\_ cil ), la temperatura de la superficie T y pared \_int el
coeficiente pelicular hc\_gas ) como de las externas (temperatura del aire exterior Text , coeficiente pelicular de transferencia de calor entre superficie exterior del cilindro y el aire exterior hc \_ext ) y las propiedades de la pared interna del cilindro (conductividad térmica Kpared, espesor e, calor específico Cp).
Se incluye, en este modelo, el efecto capacitivo que posee la pared del cilindro para el flujo de calor a través de ella.
El fenómeno de transferencia de calor a través de la pared se describe por medio de la ecuación de difusión de calor unidimensional (F. Kreith, 2001; Fernández Diéz, 2002),
ρ
Cp
∂T ∂t
= K pared
∂ 2T ∂x 2
(21)
El flujo de calor entre el gas encerrado dentro del cilindro y su pared interna es el mostrado en la ecuación (16), mientras que el flujo de calor entre la superficie externa del cilindro y el aire exterior es:
( ) •
Qcil \_ ext = Ω .h . T cil \_ ext c \_ aire cil \_ ext − Text
(22)
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
1324
G. BOURGES, J. ELIACH, M. MEDINA
•
donde Qcil \_ext representa el flujo de calor, Ωcil \_ext es la superficie externa del cilindro,
hc\_aire es el coeficiente pelicular de trasferencia de calor entre el aire y el cilindro, y
Tcil \_ ext yText representan las temperaturas de la superficie externa del cilindro y del aire exterior
respectivamente. Se discretiza la sección transversal del cilindro en una serie de secciones (Figura 8), y en
el centro de cada una cada una se sitúa un nodo ficticio en el cual se supone están concentradas las propiedades térmicas de cada sección (características del material y temperatura). Se incluyen además dos nodos exteriores a la pared, uno dentro del cilindro que representa las propiedades del gas, y el otro en el exterior del cilindro que representa las propiedades del aire exterior.
Luego, utilizando el método de las diferencias finitas (MDF) aplicado a la conducción de calor en régimen transitorio (F. Kreith, 2001; Fernández Diéz, 2002), se calculan las temperaturas de cada uno de esos nodos en función de las temperaturas del gas (Tcil ) y la
temperatura del aire exterior (Text ) con las condiciones de contorno convectivas como la dada
por las ecuaciones (16) y (22). Cabe aclarar que esta última temperatura (Text ), así como el
coeficiente pelicular hc\_aire se consideran constantes.
Gas Tcil T1 T2
T n-1 T n Text AIRE (Atmósfera)
Dx espesor de pared
Figura 8 – Pared transversal del cilindro. La misma se dividió en “n” secciones, y a cada una de estas secciones se le asignó un nodo.
Para la solución de estas temperaturas se utilizaron las ecuaciones implícitas del método de diferencias finitas debido a la mayor estabilidad que presenta dicha formulación respecto de la forma explícita. El sistema de ecuaciones queda de la forma:
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
Mecánica Computacional Vol XXV, pp. 1313-1334 (2006)
1325
T
t +Λt 1
=
2.F0 .(T
t +Λt 2
+
1 + 2.F0
Bicil .T t+Λt cil + 2.Bicil .F0
)
+T
t 1
T t+Λt 2
=
F0 .(T
t +Λt 1
+T
) t + Λt 3
+
T
t 2
1 + 2.F0
T t+Λt 3
=
F0 .(T t+Λt 2
+ T t+Λt 4 ) + T t 3
1 + 2.F0
.
.
.
.
.
.
T t+Λt n−1
=
F0 .(T t+Λt n−2 + T t+Λt n ) + T t n−1 1 + 2.F0
T t+Λt n
=
2.F0 .(T
t +Λt
n−1
+
Biext
.T
t +Λt ext
)
+
T
t n
1 + 2.F0 + 2.Biext .F0
(23)
Donde F0 es el número de Fourier, dado por:
F0
=
K pared .Δt ρ.C pared .(Δx)2
(24)
BiB cil y BBiext los números de Biot para la superficie interna y la externa del cilindro, dados por:
Bicil
=
hc \_ gas.Δx k
Biext
=
hc \_ ext .Δx k
(25) (26)
Renombrando A1 =1 + 2.F0 + 2.F0.Bicil , 2.F0
reordenando (23), queda del siguiente modo:
A2 =1 + 2.F0 + 2.F0.Biext , 2.F0
B =1 + 2.F0 , y F0
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
1326
G. BOURGES, J. ELIACH, M. MEDINA
A1.T t+Λt1 − T t+Λt 2 = Bicil .T t+Λt cil + A1.T t1 B.T t+Λt 2 − T t+Λt1 − T t+Λt 3 = B.T t 2 B.T t+Λt 3 − T t+Λt 2 − T t+Λt 4 = B.T t 3
.
.
.
.
.
.
B.T t+Λt n−1 − T t+Λt n−2 − T t+Λt n = B.T t n−1
A2.T
t +Λt n
−T
t + Λt n−1
=
B .T t+Λt
iext
extl
+
A2.T
t n
y en forma matricial:
(27)
⎡ A1 −1 0 . . 0 0 ⎤
⎢⎢−1 B −1 .
.
0
0
⎥ ⎥
⎢ 0 −1 B −1 . 0 0 ⎥
⎢ ⎢
.
...
.
⎥ ⎥
⎢.
. . . .⎥
⎢
⎥
⎢ 0 . . . −1 B −1⎥
⎣⎢104404404420 440 4−4144A32⎥⎦
M
×
⎡T
t
+Δt 1
⎤
⎢ ⎢T
t
+Δt
2
⎥ ⎥
⎢⎢T t+Δt 3
⎥ ⎥
⎢.
⎥
⎢⎢.
⎥ ⎥
⎢⎢T
t
+Δt
n−1
⎥ ⎥
⎢⎣1T4t+2Δt n43⎦⎥
T
=
⎡ ⎢ ⎢
Bicil .T
t + Δt cil
B.T
+
t 2
A1.T
t n
⎤ ⎥ ⎥
⎢ ⎢
B.T t 3
⎥ ⎥
⎢
.
⎥
⎢ ⎢
.
⎥ ⎥
⎢
B.T t n−1
⎥
⎢⎣1B4iext4.T t4+4Δt2ext 4+ A424.T4t3n ⎥⎦
N
(28)
Es decir, de la forma: M (n, n) ×T (n) = N(n) . Esto se resuelve haciendo:
T (n) = M −1(n, n) × N (n)
(29)
De esta forma se obtiene el vector T(n) cuyos componentes son las temperaturas de cada uno de los nodos del modelo. La temperatura de la superficie interna del cilindro Tpared \_int
corresponde a la primer componente del vector de temperaturas T en (28). Es decir, Tpared \_ int = T (1) .
•
2.3 Cálculo de Qaros .Calor generado por la fricción entre los aros de pistón y la pared
interna del cilindro.
En el siguiente ítem se plantea un modelo de disipación de energía por efecto de la fricción entre los aros del pistón y la pared interna del cilindro. El objetivo aquí planteado es calcular la energía por unidad de tiempo disipada durante el desplazamiento del pistón, la cual se supone será transferida al gas en forma de calor (Figura 9).
El modelo se basa en el concepto de funcionamiento de los aros de pistón que se utilizan tanto en compresores a émbolo como en motores de combustión interna (Heywood, 1988).
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
Mecánica Computacional Vol XXV, pp. 1313-1334 (2006)
1327
Qaro
Qaro
Pared
Spistón
Faro
Pcil
e
Pared
Pcil
Faro
FR
FR
d D
Figura 9: Corte donde se muestran el pistón con el aro de cierre, las paredes del cilindro, y el calor disipado por los aros de pistón.
En el modelo se propuso que la energía por unidad de tiempo disipada es proporcional a la fuerza friccional entre los aros y la superficie interna del cilindro (Faro) y al valor absoluto de la velocidad de avance del pistón ( S pistón ) (Ham, 1958).
Si llamamos con d al diámetro interno del aro (Figura 9), D el diámetro interno del cilindro que coincide con el diámetro externo del aro y e el espesor del aro, entonces la superficie expuesta a la presión interna del cilindro es (Figura 9):
Ωint\_ aro = π .d.e
(30)
Dicha superficie coincide, para el caso analizado, con la superficie de contacto entre el aro y la superficie interna del cilindro.
La fuerza no conservativa debida a la fricción entre el aro y la superficie interna del cilindro:
FR = Faro . faro
(31)
siendo Faro = Pcil.Ωint\_ aro y faro es el coeficiente de rozamiento entre el aro y la superficie.
En el cálculo de la fuerza normal sobre el aro Faro se tuvo en cuenta solo la presión interna
del cilindro. Si bien esta última presión actúa sobre el aro durante el proceso de compresión,
mientras que durante la expansión actúa la presión reinante en la parte inferior del cilindro
(menor a la existente en el cilindro), este efecto fue considerado despreciable y se consideró
solo la presión interna del cilindro.
Luego, la energía disipada por cada aro por unidad de tiempo es:
.
Q aro = −Faro . S pistón = −Pcil .Ωint\_ aro . faro . S pistón
(32)
El término de velocidad de pistón en el miembro derecho de la igualdad (32) se encuentra en valor absoluto ( S pistón ) pues el calor disipado no depende del sentido del sentido de avance del pistón, sino de su magnitud.
Entonces, para un número de aros de pistón naros , la potencia por unidad de tiempo total disipada en la fricción entre los aros y la superficie interna del cilindro es:
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
1328
G. BOURGES, J. ELIACH, M. MEDINA
.
Qaros = −Ωint\_ aro . faro .naros .Pcil . S pistón
(33)
2.4 Cálculo de ΔH .Modelo de flujo másico y entálpico.
En este ítem se plantea como objetivo obtener un modelo de flujo másico de gas a través de las válvulas de admisión y escape. Dicho flujo ocurre durante los períodos en los que cada una de ellas permanece abierta. Del balance térmico (3) se observa la necesidad de conocer el caudal de gas entrante ó saliente del cilindro, tanto para el cálculo de la variación de energía interna de la masa de gas acumulada en el cilindro, como para la variación de entalpía debida a dicho flujo másico.
El modelo utilizado para describir el flujo es basado en ecuaciones de flujo compresible a través de una restricción (Heywood, 1988; Nigro, 1999). Este modelo proviene de realizar un análisis isentrópico unidimensional, corrigiendo los efectos de gases reales mediante coeficiente un coeficiente de descarga CD.
El flujo másico se relaciona con las condiciones de estancamiento aguas arriba (p0, T0) y la presión estática aguas debajo de la restricción (la válvula), que se asume igual a la presión en la restricción pT y el área característica de la válvula Aválvula (Heywood, 1988; Nigro, 1999).
El esquema para calcular el caudal de masa que ingresa al cilindro (válvula de admisión abierta) ó sale del mismo (válvula de escape) es el siguiente:
V.A.
V.E.
P.M.S.
P.M.I.
Figura 10. Flujo a través de las válvulas.
Para flujo subsónico,
k
mientras
pT p0
>
⎡ ⎣⎢γ
2 ⎤ k −1 + 1⎦⎥ el flujo másico es
•
m
=
CD .AValvula. p0 R.T0
.(
pT p0
1
)k
.⎧⎪⎨ ⎪⎩
k2−.k1.⎢⎡⎣1
−
(
pT p0
(
)
k
−1) k
⎤⎫⎪1 ⎥⎬
/
2
⎦⎪⎭
(34)
En caso contrario, se produce choque sónico, siendo el flujo másico de la forma:
k +1
•
m
=
CD.AValvula. p0
.
R.T0
k
.⎨⎧⎩
k
2 +
⎫ 1⎭⎬
2.(
k
−1)
(35)
La variación de entalpía es entonces:
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
Mecánica Computacional Vol XXV, pp. 1313-1334 (2006)
1329
•
ΔH = m .C p .T0
(36)
Luego, la masa total encerrada en el cilindro en cada paso de tiempo es calculada como una integración a lo largo de cada ciclo.
2.4.1 Modelo de movimiento de válvulas de admisión y escape
El objetivo de este ítem es calcular la sección de paso del flujo de gas a través de las válvulas de admisión y de escape para resolver las ecuaciones de flujo másico (34) y (35) y de flujo entálpico (36). Para ello, se plantean los modelos de movimiento de las mismas .
Para simular la apertura y cierre de las válvulas de admisión y escape, se adopta un modelo de válvula de tipo masa resorte amortiguador, como el mostrado en la Figura 11.
Pexterior
V.A.
ΔP*Ωválvula
Ω
Fresorte
V.E.
Pexterior Fresorte
Ω
Pcil
Pcil
ΔP*Ωválvula
(a)
(b)
Figura 11 – (a) Esquema de válvula de admisión. (b) Esquema de válvula de escape.
Este modelo se basa en la hipótesis de que las válvulas (de aspiración y escape) funcionan por diferencia de presión entre el exterior (tubería de aspiración ó tubería de escape) y la presión interna del cilindro. Es decir, cuando la presión dentro del cilindro sea inferior a la presión de admisión ( Padm ) y la fuerza resultante de la diferencia entre dichas presiones
aplicadas a la válvula de admisión sea mayor a la fuerza que le ejerce el resorte de válvula, se abre la válvula de admisión. La misma se mantiene abierta mientras que la diferencia de presiones sea suficiente como para vencer la fuerza que ejerce el resorte de válvula sobre ella. Este proceso ocurre mientras el pistón se desplaza entre el P.M.S. (Punto Muerto Superior ) y el P.M.I. (Punto Muerto Inferior), es decir durante la carrera de descompresión o expansión.
Cuando la presión dentro del cilindro supera la presión de escape ( Pescape ) y la fuerza
resultante de la diferencia entre dichas presiones aplicadas a la válvula de escape sea mayor que la fuerza que ejercida por el resorte de la válvula de escape, se abre esta última. Esto ocurre mientras la diferencia de presiones sea superior a la fuerza ejercida por el resorte de válvula, lo cual se da mientras el pistón se desplaza entre el P.M.I. (Punto Muerto Inferior ) y el P.M.S. (Punto Muerto Superior), es decir durante la carrera de compresión. La ecuación de movimiento de las válvulas (admisión y escape) es (Ham, 1958; Edwards, 1993):
⎛⎜ ⎝
M
valvula
••
x
⎟⎞
⎠
+
⎛⎜ ⎝
R
•
x
⎟⎞ ⎠
+
(Kx)
=
F
(t)
(37)
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
1330
G. BOURGES, J. ELIACH, M. MEDINA
Donde x = x (t) es el posición de la válvula respecto a un sistema de referencia ubicado en la tapa de válvulas, y K, R y M son la rigidez del resorte de válvula (de admisión o escape), la constante de amortiguación de la válvula (que depende del tipo de válvula considerada) y la masa en movimiento de la válvula respectivamente. El término del F(t) es la excitación del sistema, y para este modelo esta dada por la diferencia de presiones entre el gas encerrado en el cilindro y la tubería de escape ó admisión (según como sea esa diferencia de presiones), multiplicada por la superficie de válvula sometida a dicha diferencia ( Ωválvula ). Es
decir,Cuando Pcil > Pesc entonces se evaluará la diferencia de presiones entre el gas dentro del cilindro (Pcil) y la existente en el conducto de escape (Pesc).
Cuando Pcil < Padm entonces se evaluará la diferencia de presiones entre el gas dentro del cilindro (Pcil) y la existente en el conducto de admisión (Padm).
Llamando con Pexterior a la presión en el conducto de escape ó admisión según sea la condición dentro del cilindro, puede escribirse la diferencia de presiones como:
ΔP = Pcil − Pexterior
(38)
Esta diferencia de presiones dependerá de la presión dentro del cilindro, la cual es una función del tiempo t, ya que la presión Pexterior se considera constante durante el período de funcionamiento de cada válvula. Por lo tanto (38) queda:
ΔP(t) = Pcil (t) − Pexterior
(39)
La presión Pcil(t) es obtenida en forma numérica por el simulador en cada paso de tiempo.
3 RESULTADOS NUMÉRICOS
En esta sección se muestra la utilización del simulador completo con condiciones de diseño impuestas como se observan en las Tabla 1 y Tabla 2.
En las figuras se muestran los resultados obtenidos luego de realizar una corrida de 10 ciclos de compresión a un régimen de 1000 C.P.M. (Ciclos Por Minuto). En la Figura 12 se ve el ciclo de compresión, diagrama P-V, en el cual se observa el pico de presión al final de la etapa de compresión politrópica y comienzo de la apertura de válvula. En las Figura 13, 14 y 15 se observan la temperatura, presión y densidad en función del tiempo que tarda el compresor en realizar cada ciclo.
Característica Caudal de entrada Presión de entrada Presión de salida Temperatura de entrada Gas Constante adiabática (γ) Constante del metano (R) Cp Cv
Valor 0,551 7x105 23x105 293 Metano (CH4) 1.321
518,31 2200 1682
Tabla 1. Características de gas.
Unidad m3/min N/m2 N/m2 ºK
J/(kg*K) J/(kg*K) J/(kg*K)
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
Mecánica Computacional Vol XXV, pp. 1313-1334 (2006)
1331
Característica Etapa nº Carrera de compresión Diámetro cilindro Espesor de pared Nº de nodos de pared Material cilindro Conduct. térmica de pared (Kpared) Calor específico de pared (Cpared) Densidad de pared Nº de ciclos por minuto Nº de ciclos por corrida Rigidez resor. válv. adm. Rigidez resor. válv. esc. Amortig. válv. adm. Amortig. válv. esc. Masa válv. Admisión Masa válv. escape
Valor 1 0.03 0.153 0.005 8 Fundición de hierro 54
480
7200 1000 10 1000 5000 0 0 0.1 0.1
Tabla 2. Cilindro de compresión.
Unidad
m m m
W/(m*K)
J/(kg*K)
Kg/m3 C.P.M. ciclos N/m N/m N*seg/m N*seg/m Kg Kg
P[Mpa]
26
24
22
20
18
16
14
12
10
8
6
0
1
2
3
4
5
6
x 10-4
V[m3]
Figura 12: Diagrama Presión vs. Volumen
T[ºK]
420
Temperatura
400
380
360
340
320
300
280
260
240 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Tiempo [seg]
Figura 13: Diagrama Temperatura vs. Tiempo
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
1332
G. BOURGES, J. ELIACH, M. MEDINA
P[Mpa]
26
Presión
24
22
20
18
16
14
12
10
8
6
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Tiempo [seg]
Figura 14: Diagrama Presión vs. Tiempo
ρ[kg/m3]
13
Densidad
12
11
10
9
8
7
6
5
4
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Tiempo [seg]
Figura 15: Diagrama Densidad vs. Tiempo
4 CONCLUSIONES
En el presente trabajo se ha modelizado el comportamiento del gas durante la compresión en un compresor obteniéndose un simulador por componentes del proceso. Se ha mostrado mediante un ejemplo numérico la aplicación del simulador para la obtención de las características del gas.
Este simulador puede resultar muy útil durante el diseño de las partes del compresor. Permite la experimentación numérica con diferentes configuraciones para evaluar la incidencia de cada componente y conseguir la combinación óptima de los mismos. En particular se puede aplicar al diseño de las válvulas y de las aletas de enfriamiento. En estas últimas se utiliza el simulador para estimar la temperatura superficial del cilindro de compresión.
En futuros trabajos se puede extender el desarrollo del simulador al caso de un compresor multi-etapa y al estudio multiparamétrico del rendimiento volumétrico.
5 NOMENCLATURA
Símbolo Cantidad
t
Tiempo.
Δt
Paso temporal.
φ
Angulo girado por el cigüeñal.
m
Masa de gas.
•
m
Flujo másico de gas por unidad de tiempo.
s
Distancia entre eje de cigüeñal y perno de
pistón.
•
s
Tasa de variación de la distancia s respecto
del tiempo.
Unidades [S.I.]
seg. seg. Radianes Kg Kg/seg.
m
m/s
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
Mecánica Computacional Vol XXV, pp. 1313-1334 (2006)
1333
S pistón S pistón
a L D Vc Vcil
•
Vcil
Ωint\_ aro
AValvula
CD Pcil
E int
•
Eint
Velocidad de desplazamiento del pistón.
Velocidad promedio de desplazamiento
del pistón.
Radio del cigüeñal.
Longitud de biela.
Diámetro del cilindro.
m
Volumen nocivo encerrado en el cilindro.
Volumen encerrado en el cilindro.
Variación del volumen encerrado dentro del cilindro respecto del tiempo. Superficie del aro en contacto con la presión
del cilindro. Sección de paso del flujo a través de cada válvula. Coeficiente de descarga de la válvula. Presión de la masa de gas encerrada en el cilindro. Energía interna.
Variación de energía interna en el tiempo.
W
•
W
•
Q pared
•
Q aros
ΔH f aro
Kgas Kpared k R Cp Cv Cpared Re Pr BiB F0 hc
Trabajo transferido.
Potencia; trabajo por unidad de tiempo.
Flujo de calor por unidad de tiempo
a través de la pared del cilindro.
Flujo de calor por unidad de tiempo
debido al rozamiento entre los aros de pistón y la superficie interna del cilindro. Flujo entálpico a través de las válvulas. Coeficiente de fricción entre cada aro y la
superficie del cilindro. Coef. de conductividad térmica del gas. Coef. de conductividad térmica de la pared. Coeficiente politrópico del gas. Constante del gas. Calor específico a presión constante del gas. Calor específico a volumen constante del gas. Calor específico del material del cilindro. Número de Reynolds. Número de Prandtl. Número de Biot. Número de Fourier. Coeficiente pelicular de transferencia de calor.
m/s m/s
m m
m m3 m3/seg.
m2
m2
--N/m2
J
J/seg.
J J/seg. J/seg.
J/seg.
---
J/mºK J/mºK ----J/kgºK J/kgºK J/kgºK --------J/m2ºK
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
1334
G. BOURGES, J. ELIACH, M. MEDINA
REFERENCIAS
R. B. Bird, W. E. Steward, E. N. Lightfood. Fenómenos de transporte. Editorial Reverté, S.A., ISBN 84-291-7050-2, 2001.
C. H. Edwards, Jr.; D. E.Penney. Ecuaciones diferenciales elementales y problemas con condiciones en la frontera. Tercera edición. Prentice-Hall Hispanoamericana, S.A, ISBN 968-880-414-2. 1993.
V.M. Faires. Termodinámica. Editorial LIMUSA, ISBN 968-18-3943-9, 1999. P. Fernández Diéz. Ingeniería Energética. Departamento de Ingeniería Eléctrica y
Energética. Universidad de Cantabria, 2002. C. W. Ham, E. J. Crane, W. L. Rogers. Mechanics of machinery. Fourth Edition.
International Student Edition, McGraw-Hill, INC. 1958. J. B. Heywood. Internal Combustion Engine Fundamentals. McGraw-Hill. Series in
Mechanical Engineering, ISBN 0-07-028637-X, 1988. J.B. Jones and R.E. Dugan. Ingeniería Termodinámica. Primera edición. Prentice-Hall
Hispanoamericana, S.A, ISBN 968-880-845-8, 1997. F. Kreith, M. S. Bohn. Principios de la transferencia de calor. Sexta edición. Thomson
Learning, ISBN 970-686-063-0 , 2001. S. Nakamura. Analisis numérico y visualización gráfica con MATLAB®. Primera
edición. Prentice-Hall Hispanoamericana, S.A, ISBN 968-880-860-1, 1997. N. Nigro, M. A. Storti, L. Ambroggi. Modelización numérica de un motor de
combustión monocilíndrico encendido por chispa. Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería. ISSN 0213-1315. Volumen 15, 1, 21-54, 1999. L.F. Shampine and M.W. Reichelt. The MATLAB ODE suite. On line Manuals of
MATLAB®.
Copyright © 2006 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar
Ver+/-
