Embed Size (px)
Citation preview
UNIVERSIDADE FEDERAL DO RIO DE JANEIRO Departamento de Engenharia Mecânica
DEM/POLI/UFRJ
SIMULAÇÃO COMPUTACIONAL DA COMBUSTÃO EM UM MOTOR DIESEL
MARÍTIMO DE ROTAÇÃO CONSTANTE
Filipe Augusto Serrão Matias
PROJETO FINAL SUBMETIDO AO CORPO DOCENTE DO DEPARTAMENTO DE
ENGENHARIA MECÂNICA DA ESCOLA POLITÉCNICA DA UNIVERSIDADE
FEDERAL DO RIO DE JANEIRO COMO PARTE DOS REQUISITOS
NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE ENGENHEIRO MECÂNICO.
Aprovado por:
________________________________________________
Prof. Marcelo José Colaço
________________________________________________
Prof. Albino José Kalab Leiróz
________________________________________________
Prof. Hélcio Rangel Barreto Orlande
________________________________________________
Prof. Nome do Membro da Banca Opcional
RIO DE JANEIRO, RJ - BRASIL
JULHO DE 2014
ii
AGRADECIMENTOS
Ao meu orientador Prof. Marcelo José Colaço pela oportunidade que me foi dada
em realizar este projeto e por todo suporte dado durante o desenvolvimento.
Aos professores Albino José Kalab Leiróz e Hélcio Rangel Barreto Orlande, pela
disponibilidade em participarem da minha banca examinadora.
À minha família, principalmente aos meus pais, por todo apoio e suporte que
sempre me foi dado.
À minha namorada Renata Cardoso Bezerra pelo carinho e apoio nos momentos
mais difíceis.
À meu grande amigo Yuri Antelo Barbosa, por todas as vitórias que
conquistamos nestes últimos anos.
iii
RESUMO
Este projeto tem por objetivo a simulação computacional de um motor tipo
diesel, para estudo da variação do consumo de combustível de diferentes misturas de
diesel, biodiesel e etanol. A partir de dados de pressão, de ensaios experimentais do motor
com diesel marítimo, os parâmetros da simulação foram ajustados. Utilizou-se os
modelos de Wiebe duas zonas para simular a combustão, e o modelo desenvolvido pela
empresa AVL em 2000 para estimar o coeficiente de transferência de calor entre os gases
de combustão e as paredes do cilindro. A partir das curvas de pressão, foram calculados
o IMEP, a pressão máxima e seu respectivo ângulo, além de ser proposta uma
metodologia para calcular o início e o fim da combustão.
iv
ABSTRACT
This Undergraduate Project presents a diesel engine combustion model, in which
different blends of ethanol, diesel and biodiesel were simulated, and their respective fuel
consumption analysed. Using experimental data of the in cylinder pressure for diesel, the
parameters were adjusted and the model validated. The Wiebe two zones combustion
model was used to simulate heat release; and to simulate the heat transfer coefficient, a
model created by AVL and published in 2000 was used. From the pressure experimental
data, the IMEP, and peak pressure were calculated; also, it was proposed a method to
calculate the combustion start and its duration.
v
Sumário
1. Introdução ............................................................................................................................. 1
2. Revisão bibliográfica............................................................................................................. 2
2.1. Misturas de óleo diesel e biodiesel ................................................................................ 2
2.2. Misturas de óleo diesel, biodiesel e etanol .................................................................... 5
2.3. Simulação de Combustão .............................................................................................. 6
2.4. Modelos de Combustão ................................................................................................. 7
3. Motores de combustão interna .............................................................................................. 9
3.1. Funcionamento dos motores à combustão interna ........................................................ 9
3.2. Classificação de motores ............................................................................................. 12
4. Metodologia ........................................................................................................................ 16
4.1. Motor MAN Innovator 4C .......................................................................................... 16
4.2. Combustíveis ............................................................................................................... 17
4.2.1. Diesel ................................................................................................................... 20
4.2.2. Etanol .................................................................................................................. 22
4.2.3. Biodiesel .............................................................................................................. 26
4.2.4. Misturas de combustíveis .................................................................................... 29
4.3. Modelo de combustão ................................................................................................. 29
4.3.1. Modelo de combustão de Wiebe ......................................................................... 32
4.4. Modelo de transferência de calor ................................................................................ 37
4.5. Programa para análise das curvas de pressão experimentais ....................................... 39
4.5.1. Início da combustão ............................................................................................ 40
4.5.2. Final da combustão .............................................................................................. 46
4.5.2.1. Método do erro ............................................................................................ 46
4.5.2.2. Método da menor derivada de ln𝑃 .............................................................. 48
4.5.3. Cálculo de IMEP ................................................................................................. 50
4.5.4. Cálculo de BMEP e FMEP .................................................................................. 52
4.5.5. 𝑃𝑚𝑎𝑥 e 𝜃𝑃𝑚𝑎𝑥 .................................................................................................. 53
4.5.6. Tratamento dos dados .......................................................................................... 53
4.5.7. Curva de pressão representativa do ciclo ............................................................ 54
4.5.8. Resultados do programa das curvas de pressão ................................................... 55
4.6. Cálculo da posição das válvulas de admissão e exaustão............................................ 57
5. Construção do modelo no BOOST ...................................................................................... 63
5.1. Ajustes ......................................................................................................................... 64
5.1.1. Turbo compressor ................................................................................................ 64
vi
5.1.2. IMEP ................................................................................................................... 65
5.1.3. FMEP .................................................................................................................. 66
6. Resultados ........................................................................................................................... 68
6.1. Validação da simulação ............................................................................................... 68
6.2. Consumo de combustível para diferentes misturas ..................................................... 71
6.2.1. 100% de Carga nominal do motor ....................................................................... 71
6.2.2. 50% de carga ....................................................................................................... 73
7. Conclusões .......................................................................................................................... 75
8. Referências Bibliográficas ................................................................................................... 76
Anexo I: Programa para calcular dados da curva de pressão e utilizando o método descrito no
item 4.5.2.1. ................................................................................................................................. 82
Anexo II: Programa para calcular dados da curva de pressão e utilizando o método descrito no
item 4.5.2.2. ............................................................................................................................... 103
Anexo III: Calculo da abertura das válvulas de admissão e exaustão ....................................... 121
Anexo IV: Programa para calcular o deslocamento das válvulas ............................................. 122
1
1. Introdução
Na Europa, a nova legislação EURO VI já entrou em vigor, com limites de
emissões de poluentes mais restritos que os determinados pela EURO V que a antecedeu.
No Brasil a nova fase do CONAMA, a P8, deve entrar em vigor até 2016 e acredita-se
que trará cortes de emissões similares aos da legislação EURO VI [1].
A crescente preocupação mundial com relação à poluição atmosférica e as novas
legislações, como as citadas, que surgem em paralelo, são alguns dos fatores responsáveis
pelos altos investimentos em pesquisa no setor automobilístico.
Além de soluções imediatas, países como os da União Europeia, os quais vêm
se planejando a longo prazo, assinaram um acordo em 2007 se comprometendo a que
20% do consumo de energia seja de fontes renováveis [2]. O acordo também prevê que
pelo menos 10% da energia utilizada no setor de transporte seja proveniente de fontes
renováveis. Especialistas esperam que boa parte da cota do setor automotivo seja suprida
pelo aumento do uso de biocombustíveis [3]
O uso de biocombustíveis em larga escala já é uma realidade em diversos países,
e a tendência é um aumento significativo no seu consumo [4]. Além de substituir em parte
o consumo de combustíveis fósseis, o uso de biocombustíveis traz outras vantagens, como
é o caso do etanol, responsável pelo aumento da octanagem da gasolina [5], e do biodiesel
que tende a diminuir as emissões de MP (material particulado), CO (monóxido de
carbono) e HC (hidrocarbonetos) [6].
A modelagem de um motor a combustão interna é de extrema complexidade, e
ainda hoje não é totalmente dominada. Como ela é altamente vantajosa por apresentar
baixos custos e curtos tempos, em comparação com os testes experimentais, inúmeras
pesquisas estão sendo realizadas nesta área, a fim de buscar simulações mais acuradas.
Entretanto, apesar dos constantes avanços encontrados, muito ainda precisa ser feito,
principalmente nas áreas de simulações com diferentes composições dos combustíveis.
Há diversas pesquisas envolvendo combustíveis com diesel misturado a etanol,
como é o caso das realizadas por ROCHA [7], SHI et al. [8] e AL-HASSAN et al. [9].
Porém, não foi encontrado estudos de simulação da combustão em motores envolvendo
esses tipos de misturas. Sendo assim, foi escolhido como assunto deste trabalho a análise
computacional do consumo de diferentes misturas de diesel, biodiesel e etanol.
2
2. Revisão bibliográfica
2.1. Misturas de óleo diesel e biodiesel
A mistura de óleo diesel com biodiesel já é utilizada em muitos países do mundo
[4]. Esta mistura apresenta muitas vantagens, como diminuição do consumo de
combustíveis provenientes do petróleo, decréscimo das emissões de gases responsáveis
pelo efeito estufa [10] e redução das emissões de certos poluentes veiculares [6]. Muitos
trabalhos já foram realizados com misturas desses dois combustíveis e alguns serão
comentados a seguir.
KARAVALAKIS et al. [11] realizaram experimentos com biodiesel proveniente de
diferentes matérias primas a fim de identificar alguma correlação entre as emissões e o
tipo de biodiesel. Misturas de diesel com percentuais de 5%, 10% e 20% de biodiesel
foram testadas, sendo que os biodieseis utilizados foram os seguintes: éster metílico a
base de soja; éster metílico a base de palma misturada com óleo de côco; éster etílico a
base de colza misturado com óleo de girassol e resíduos de óleo de cozinha.
O veículo utilizado para os ensaios foi o Hyundai i-10 (CRDi VGT) com tecnologia
common-rail. O ciclo Common Artemis Driving Cycles (CADC) foi utilizado para
medição de material particulado (MP), monóxido de carbono (CO), dióxido de carbono
(CO2), Óxidos de Nitrogênio (NOx) e hidrocarbonetos (HC). Já o ciclo New European
Driving Cycle (NEDC) foi utilizado para medir as emissões de MP, NOx, e CO2.
Os autores não encontraram uma correlação entre as emissões de CO2 e os
percentual de biodiesel usado. Enquanto as emissões de CO, HC e MP diminuíram com
o aumento de biodiesel na mistura, as de NOx aumentaram. Além disso não foi
encontrada correlação entre as emissões e o tipo de biodiesel.
Em 2002 o EPA (Environmental Protection Agency) publicou um relatório técnico
[6] no qual foram analisados mais de 300 dados de trabalhos distintos com o intuito de se
ter uma visão geral dos impactos do uso do biodiesel nas emissões veiculares. Entre os
diversos estudos realizados, um deles foi a criação de uma correlação entre a quantidade
de biodiesel no combustível e a mudança percentual das emissões. A correlação utilizou
a equação 1 para descrever as emissões de NOx, MP, HC e CO.
3
% 𝑑𝑒 𝑚𝑢𝑑𝑎𝑛ç𝑎 𝑑𝑎𝑠 𝑒𝑚𝑖𝑠𝑠õ𝑒𝑠
= {−1 + exp[𝑎 ∗ (𝑣𝑜𝑙% 𝑑𝑒 𝑏𝑖𝑜𝑑𝑖𝑒𝑠𝑒𝑙)]} ∗ 100%
(1)
Os valores da constante 𝑎 foram ajustados utilizando os dados experimentais. A
tabela 1 apresenta tais valores.
Tabela 1: Valores da constante "a" para a equação 1
Constante 𝒂
NOx 0,0009794
MP -0,006348
HC -0,011195
CO -0,006561
Substituindo os valores da tabela 1 na equação 1 tem-se um gráfico (figura 1) das
tendências das emissões de poluentes em função do percentual de biodiesel misturado ao
óleo diesel.
Figura 1: Gráfico representando as tendências das emissões de NOx, PM, CO e HC em
função do teor de biodiesel
4
ROBBINS et al. [12] analisaram diversos trabalhos já publicados com relação às
emissões de misturas de diesel e biodiesel. Estes trabalhos realizaram experimentos em
diversos motores e veículos pesados de diferentes fabricantes e nacionalidades.
Os dados provenientes de experimentos em motores pesados permitiram que os
autores concluíssem que há tendência de diminuição das emissões de HC, CO e MP
conforme aumenta-se a quantidade de biodiesel na mistura. Pode-se também concluir que
as emissões de NOx tendem a aumentar com o aumento do teor de biodiesel, apesar desse
aumento ter uma taxa de crescimento menor que a taxa de diminuição das emissões dos
outros poluentes. Foram também analisadas as emissões de veículos pesados e, apesar de
não serem iguais, os resultados são similares aos testes em motores.
Os autores realizaram uma comparação entre os dados encontrados por eles e o
gráfico publicado pelo EPA em 2002 [6]. A figura 2 apresenta um gráfico comparando as
emissões de NOx, CO, HC e MP em função do percentual de biodiesel misturado ao óleo
diesel. As linhas pontilhadas representam a tendência dos dados encontrados por
ROBBINS et al. [12] enquanto as linhas contínuas representam o gráfico original
publicado pela EPA.
Figura 2: Comparação entre os dados de emissões do EPA [6] e de ROBBINS et al.
[12]
5
2.2. Misturas de óleo diesel, biodiesel e etanol
A mistura ternária de óleo diesel, biodiesel e etanol já é estudada há algum
tempo. Duas das grandes vantagens do uso de biocombustíveis são a diminuição do
consumo de combustíveis fósseis e a diminuição das emissões de gases responsáveis pelo
efeito estufa. Apesar dessas vantagens, o grande interesse no estudo de misturas com
etanol em motores diesel é a tendência de diminuição das emissões de NOx. A seguir
serão citados alguns trabalhos que abordaram essas misturas.
ACYOLY et al. [13] testaram formulações de óleo diesel com adição de
biodiesel e etanol em diversas proporções a fim de verificar suas propriedades físico-
químicas. Posteriormente as misturas foram qualificadas quanto à possibilidade de
utilização em motores diesel.
Uma das vantagens em se utilizar misturas de diesel com biodiesel para adicionar
etanol é que o biodiesel tem afinidade química tanto com o diesel quanto com o etanol,
facilitando a solubilidade da mistura. Apesar da afinidade do biodiesel com o etanol, para
as misturas utilizadas pelo autor (10%, 15% e 20% de etanol), o biodiesel não foi
suficiente para manter a estabilidade da mistura. Portanto, um tensoativo foi utilizado
para manter a mistura homogenia.
Todas as combinações de diesel com 5, 10, 15 e 20% de biodiesel; 10, 15 e 20%
de etanol; e 1% do tensoativo foram criadas e então analisadas. Com exceção da mistura
de 20% de biodiesel e 20% de etanol, todas as outras se enquadraram na resolução n° 7
da ANP de 2008 [14] no que diz respeito ao índice de cetano, viscosidade, massa
específica e estabilidade.
Os autores identificaram que na ausência do tensoativo, para teores de 5% de
biodiesel, foi possível a solubilidade de até 1% de etanol; e para teores de 10% de
biodiesel, solubilidades pouco maiores que 2%.
AL-HASSAN et al. [9] analisaram a estabilidade e o desempenho de misturas
ternárias (diesel, biodiesel e etanol anidro). Mantendo fixa a concentração de biodiesel
em 10%, misturas com 5, 10, 15 e 20% de etanol foram realizadas. A mistura com 5% de
etanol foi a que apresentou melhor estabilidade. Utilizando um motor de 6 kW de
potência, a rotação foi variada entre 800 e 1600 rpm com intervalos de 200 rpm, sendo
que para cada condição do motor foram analisados o diesel puro (combustível base) e as
6
misturas realizadas. A mistura com 5% de etanol apresentou resultados bastante
semelhantes ao combustível base, tanto em potência quanto em consumo de combustível.
Notou-se que conforme era aumentado o teor de etanol, a potência produzida diminuía e
o consumo aumentava. Concluiu-se que pequenos teores de etanol podem ser utilizados
em motores diesel sem que haja necessidade de mudanças no projeto. Os autores
enfatizaram que dentre os combustíveis testados, o que apresentou os melhores resultados
à parte do diesel, em termo de eficiência do motor, foi o de 5% de etanol e 10% de
biodiesel.
ROCHA [7] realizou um estudo de misturas ternárias de combustíveis (diesel,
biodiesel e etanol). Dentre outras análises, foram verificados, para as misturas que se
mantiveram estáveis, o número de cetano e o atraso de ignição em um motor CFR/ASTM-
cetano. Diante dos combustíveis analisados encontrou-se que, para qualquer mistura
contendo apenas diesel e biodiesel, e para misturas com até 5% de etanol com no mínimo
10% de biodiesel, o número de cetano atende à resolução n° 50 de 2013 da ANP [15].
Misturas com até 5% de etanol e até 100% de biodiesel, quando submetidas às mesmas
condições de operação do motor, variaram menos de 1° no atraso de ignição quando
comparadas com o combustível base (5% de biodiesel e 95% de óleo diesel). Também
observou-se, que misturas contendo etanol apresentaram redução das emissões de NOx
em comparação ao combustível base.
2.3. Simulação de Combustão
Apesar de nos últimos anos muitos estudos terem sido realizados na área de
simulação de combustão em motores de combustão interna, ainda existem muitos desafios
a serem vencidos. A área de simulação de combustão com ênfase nas propriedades dos
combustíveis é uma das áreas deficitárias, precisando ainda de maiores estudos [5]. A
seguir serão apresentados alguns trabalhos realizados sobre a simulação da combustão
em motores.
LOGANATHAN et al. [16] buscaram simular a combustão e a performance de
um motor diesel turbo alimentado, abastecido de uma mistura de diesel com um
combustível oxigenado. Para simular a combustão foi utilizado o modelo de combustão
de Wiebe duas zonas [17]. No caso da transferência de calor pela câmara de combustão
foi utilizado o modelo de ANNAND [18]. Dois combustíveis foram analisados nesse
7
trabalho: um diesel puro comercial e uma mistura desse combustível puro com 5% v/v de
diglicol metil éter - éter bis(2-metoxietílico) - que é o combustível oxigenado.
Os dados simulados apresentaram pequena diferença, para ambos os
combustíveis, quando comparados com dados experimentais do motor. Sendo assim o
modelo proposto pelos autores pôde ser validado. O combustível oxigenado apresentou
menor consumo de combustível e maior eficiência térmica quando comparado com o óleo
diesel comercial.
MELO [5] realizou a simulação de um motor FLEX utilizando o programa
BOOST da empresa AVL. Foram utilizadas diferentes misturas de gasolina comercial
(contendo 25% de álcool anidro) com etanol hidratado. Foi utilizando o modelo de
combustão de Wiebe duas zonas para simular a liberação de calor durante a combustão,
e o modelo de Woschni publicado em 1967 [19] para calcular o coeficiente de
transferência de calor da mistura.
Os parâmetros do modelo de Wiebe foram ajustados individualmente para cada
dado experimental, e os parâmetros de Woschni foram mantidos constantes para cada
condição do motor, independentemente do combustível.
Para a simulação das emissões houve necessidade de se ajustar os parâmetros
para cada mistura de combustível e cada condição do motor. Concluiu-se então que o
software de simulação utilizado (BOOST) não é preditivo para simulação de emissões de
poluentes.
2.4. Modelos de Combustão
Há diversos modelos de combustão para motores de combustão interna. Estes
modelos podem ser agrupados de diversas maneiras. MACHADO [20] faz uma breve
descrição da classificação proposta por STONE [21] e completa com uma quarta
classificação (Algébrica) proposta por BARROS [22]. Os quatro grupos: Algébrico, Zero
dimensional, Quase dimensional e Multidimensional serão comentados a seguir [20].
Algébrico: Utiliza um sistema de equações algébricas baseadas em
relações termodinâmicas clássicas e utiliza-se de eficiências para corrigir
os resultados teóricos calculados. Requer o menor tempo computacional.
Zero dimensional: utiliza um sistema de equações diferenciais ordinárias
e as associa a dados empíricos para ajustá-las.
8
Quase dimensional: utiliza um sistema de equações diferenciais parciais
associadas a modelos semi-empíricos, como os modelos de turbulência
usados pelo modelo de combustão Fractal [23].
Multidimensional: semelhante ao modelo Quase dimensional, porém,
além da dimensão temporal, esses modelos utilizam mais de uma
dimensão física.
MELO [24] explica melhor os modelos de combustão Zero dimensionais citando
os quatro subgrupos descritos por LANZAFAME et al. [25]. Os autores propõem que os
modelos Zero dimensionais sejam divididos em quatro categorias, conforme descritas a
seguir.
k (razão de calores específicos) é considerado constante, porém não é
considerada a troca de calor entre os gases no interior do cilindro e as
paredes da câmara de combustão;
k é considerado constante, e é computado a troca de calor entre gases e
as paredes da câmara;
k varia em função da temperatura, mas não é considerada a troca de calor
entre gases e as paredes da câmara;
k varia em função da temperatura e é considerada a troca de calor entre
gases e as paredes da câmara.
Dentre as quatro categorias mencionadas acima, a última é considerada a melhor
para simular a combustão em um motor a combustão interna [24]. Como o k varia em
função da temperatura e a troca de calor é considerada, esta categoria se aproxima mais
do fenômeno físico.
9
3. Motores de combustão interna
3.1. Funcionamento dos motores à combustão interna
Os motores de combustão interna são maquinas capazes de converter em
trabalho parte da energia química acumulada em um combustível. Podemos observar na
figura 3 [17] os principais componentes de um motor de combustão interna.
Uma imagem simplificada do cilindro e do eixo de manivela é encontrada na
figura 4. Na figura pode-se observar o Ponto Morto Superior (PMS), que é o ponto
máximo que o pistão atinge durante o processo de compressão; o Ponto Morto Inferior
(PMI), que é o ponto mínimo atingido pelo pistão durante o processo de compressão; o
volume deslocado (𝑉𝑑), que é o volume deslocado pelo pistão durante o seu curso; o
Figura 3: Desenho em corte de um motor de combustão interna
Eixo de Manivelas
(Virabrequim)
10
volume morto do cilindro (𝑉𝑐), que é o volume mínimo do cilindro (ocorre quando o pistão
se encontra no PMS); e o volume total no interior do cilindro (𝑉𝑙).
Figura 4: Vista simplificada de corte lateral do cilindro incluindo o conjunto biela e
manivela
Pode-se observar na figura 4 que o movimento vertical do pistão de ida e volta
movimenta um eixo (eixo de manivelas ou virabrequim), que realiza um movimento
cíclico. É no eixo do virabrequim que o trabalho produzido pelo motor é disponibilizado.
O pistão se movimenta entre dois pontos limites: o Ponto Morto Superior (PMS) e o Ponto
Morto Inferior (PMI), sendo este intervalo denominado Curso do Pistão. Um motor de
quatro tempos tem seu funcionamento dividido em quatro etapas: admissão, compressão,
combustão (ou expansão) e exaustão. Para realizar essas quatro etapas, completando um
ciclo do motor, o virabrequim realiza duas voltas completas.
Na admissão o pistão se movimenta no sentido do PMS para o PMI, passando a
aumentar o volume da câmara. Nessa etapa a válvula de admissão é aberta e um fluido é
admitido dentro do cilindro. Este fluido é apenas ar atmosférico, caso seja um motor
PMS
PMI 𝑉
𝑑 𝑉
𝑙
𝑉 𝑐
11
Diesel; ou ar atmosférico misturado ao combustível, no caso de motores do tipo Otto. Na
etapa seguinte, compressão, as duas válvulas se encontram fechadas e o gás no interior
do cilindro é comprimido pelo pistão, que se move no sentido de PMI para PMS. A
próxima etapa, combustão, é definida por elevadas pressões no interior do cilindro devido
ao processo de combustão dos gases. Concomitantemente, o pistão expande (se
movimenta do PMS para o PMI) e é nessa etapa que o motor produz trabalho. No caso
dos motores Diesel, combustível é injetado próximo do PMS e o início da combustão
ocorre logo em seguida; já no caso de motores Otto, o combustível já se encontra
misturado ao ar no interior do cilindro, e a combustão ocorre devido a uma centelha. Na
última etapa, exaustão, o pistão se movimenta do PMI para o PMS, no sentido de
comprimir os gases, porém, como a válvula de descarga está aberta, os gases no interior
do cilindro (produtos da combustão) são expelidos do mesmo. Deve-se notar que o motor
é uma máquina cíclica, ou seja, após a realização da etapa de exaustão o motor volta à
etapa de compressão e recomeça o ciclo.
Com exceção da etapa de expansão, onde ocorre a combustão, nas outras três
etapas o motor aspirado apenas consome trabalho. Sendo assim, para que haja um balanço
positivo de geração de trabalho, a etapa da combustão deve gerar mais trabalho que todo
trabalho consumido em um ciclo do motor. Segundo HEYWOOD [17], na etapa da
combustão gera-se aproximadamente cinco vezes mais trabalho que o motor consome em
um ciclo. No caso de motores turbocomprimidos, a pressão do gás pode ser alta suficiente
a ponto dos gases gerarem trabalho sobre o pistão na etapa de admissão. As quatro etapas
de um motor a combustão interna se encontram representadas separadamente na figura 5.
12
Figura 5: vista lateral do cilindro dividida nas quatro etapas de um motor de quatro
tempos (Admissão, Compressão, Expansão e Exaustão).
3.2. Classificação de motores
Existem diversas classificações de motores à combustão interna. Algumas delas
são apresentadas abaixo [17]:
a- Aplicação: Automobilístico, locomotiva, gerador, aviões de pequeno porte,
etc.
b- Projeto do motor: Pistões em V, em linha, radial, etc.
c- Ciclo de trabalho: Motor de quatro tempos aspirado (admite ar a pressão
atmosférica) ou turbo (admite ar a pressão maior que a atmosférica); dois
tempos aspirado ou turbo; etc.
d- Combustível: Gasolina, óleo diesel, etanol hidratado, biodiesel, gás natural,
óleo combustível bunker, etc.
13
e- Método de ignição: Ignição por centelha (geralmente para combustíveis que
se assemelham à gasolina) ou ignição por compressão (geralmente para
combustíveis que se assemelham ao óleo diesel), etc.
f- Fluido de arrefecimento: Motores resfriados à água, ar, óleo, etc.
Essas são algumas das classificações de motores. Pode-se observar apenas pelas
classificações o número variável de projetos de motores de combustão interna existentes
atualmente. Habitualmente os motores automobilísticos são classificados em duas
categorias, de acordo com seu método de ignição: motor Diesel (ignição por compressão)
ou motor Otto (ignição por centelha).
Nicolaus Otto foi o primeiro a construir, em 1876 [26], um motor de combustão
interna de quatro tempos, dando origem hoje a motores do ciclo Otto (ou motores Otto).
Este motor é amplamente utilizado pela indústria automobilística para a queima de
combustíveis como gasolina e etanol. O motor se baseia no princípio dos quatro tempos
mencionados no item 3.1, sendo que na admissão, o motor admite uma mistura quase
homogênea de ar e combustível, sendo esta mistura comprimida em seguida. Para o
começo da combustão do gás, uma centelha é liberada por um mecanismo chamado de
vela e o início da queima do gás começa neste ponto e se propaga até os limites da câmara
de combustão.
Rudolf Diesel foi o inventor do que hoje é conhecido por motor ciclo Diesel (ou
motor Diesel). Em 1897, Diesel apresentou ao mundo seu motor monocilindro de 20 HP
(14,7kW) que era capaz de converter 26,8% da energia contida no combustível em
potência [27]. Assim como o motor de Nicolaus Otto, este motor também operava como
um motor de quatro tempos, porém, na admissão, apenas ar é admitido no interior do
cilindro e somente este gás participa da etapa de compressão. Próximo ao fim da
compressão, o combustível é injetado diretamente no interior do cilindro e, devido à alta
pressão e temperatura do gás, e a presença de oxidante, a combustão é iniciada.
Hoje a maioria dos motores automobilísticos operam com motores Diesel ou
Otto e, sabe-se que para estes motores, quanto maior razão de compressão do motor (razão
entre os volumes da câmara no início e no final da etapa de compressão) maior é a sua
eficiência térmica (conversão da energia química do combustível em trabalho). Os
motores Otto trabalham com um limite de razão de compressão menor que os motores
Diesel, causado por uma limitação do combustível. Para que haja um bom funcionamento
14
do motor é necessário que haja apenas um ponto de começo da combustão, que é pela
centelha da vela. Por isso a mistura de ar e combustível deve resistir às altas pressões e
temperaturas durante a etapa de compressão sem que haja autoignição da mistura. Essa
qualidade do combustível é desejada, pois se houver mais de um ponto de início da
combustão, quando a frente de chama originada desses pontos colidirem, as pressões no
interior do cilindro alcançarão valores bem maiores que aqueles para o qual o motor foi
projetado para funcionamento normal.
Dois métodos foram padronizados internacionalmente para medir a capacidade
do combustível usado nos motores Otto de suportar as altas pressões e temperaturas no
interior do cilindro. Esses métodos indicam essa qualidade de forma numérica, conhecida
como número de octano ou octanagem do combustível. O “Research Method” descrito na
norma ASTM D-2699 [28] calcula o RON (Research Octanage number), enquanto o
“Motor Method”, encontrado na norma ASTM D-2700 [29], calcula o MON (Motor
Octanage number). Esses métodos são realizados no motor padronizado CFR F-2U.
Apesar de ser possível conseguir números de octanagem maiores do que os encontrados
hoje na gasolina comercial, o custo para produção desse combustível em larga escala
aumentaria consideravelmente o valor para o consumidor final. Sendo assim, ainda é
impraticável comercializar gasolinas com octanagem consideravelmente maiores que as
encontradas hoje no mercado. Por esse motivo, a razão de compressão do motor Otto é
limitada pelo número de octano do combustível, o que restringe consideravelmente seu
rendimento térmico.
Os motores Diesel trabalham de forma diferente dos motores Otto. Para o melhor
funcionamento do motor, deseja-se que o combustível tenha a menor capacidade de
resistir às pressões e temperaturas no interior do cilindro, ou seja, deseja-se que o
combustível comece a sua queima tão logo quanto ele é injetado no cilindro. Esta
qualidade do combustível é desejada para que o mesmo, ao ser injetado na câmara,
comece o processo de combustão antes que uma grande massa de combustível se acumule
causando uma liberação de calor alta e repentina.
Analogamente aos combustíveis usados nos motores Otto, um método
normalizado utilizado internacionalmente quantifica o combustível usado nos motores
Diesel quanto ao atraso de início da combustão (intervalo de tempo entre o momento que
o combustível é injetado no cilindro e o início da combustão). Essa qualidade do
15
combustível é chamada de número de cetano ou cetanagem do combustível (ASTM
Method D613 [30]). Como no caso do motor Diesel é desejável uma baixa resistência do
combustível à combustão, é vantajoso altas temperaturas e pressões no interior do
cilindro, podendo então aumentar ao máximo a razão de compressão do motor. Neste
caso, a razão de compressão é limitada pela capacidade dos materiais que compõem o
motor a resistirem às suas condições de operação. Sendo assim, é possível construir
motores Diesel com taxas de compressão muito mais elevadas que dos motores Otto, e é
por esse motivo que os motores Diesel atingem eficiências térmicas mais elevadas.
16
4. Metodologia
4.1. Motor MAN Innovator 4C
Neste projeto foi modelado o motor modelo Innovator 4C do fabricante MAN,
o qual é um motor Diesel e apresenta potência nominal de 500 kW. Tal motor possui
rotação constante e pode operar tanto com óleo diesel quanto com óleos combustíveis
conhecidos como Bunker. Na tabela 2 são apresentados dados técnicos do motor [31].
Tabela 2: Especificação do motor da MAN modelado neste trabalho
Cilindros 5 em linha
Diâmetro do Cilindro (𝑫) 160mm
Curso do Pistão (𝑳) 240mm
Comprimento da Biela (𝒍) 480mm
Cilindrada Total 24000 cm3
Taxa de Compressão (𝒓𝒄) 15,2:1
Regime Estacionário
Rotação 1200 RPM
Sequência de combustão 1-2-4-5-3
Potência 500 kW
Válvulas por Cilindro 4 válvulas
Diâmetro das Válvulas de Admissão 57 mm
Folga das Válvulas de Admissão 0,4 mm
Diâmetro das Válvulas de Exaustão 54 mm
Folga das Válvulas de Exaustão 0,5 mm
A figura 6 [32] é uma foto do motor instalado em um banco de provas do
Laboratório de Maquinas Térmicas, o Bunker I.
17
4.2. Combustíveis
Foram utilizados três combustíveis bases para fazer as misturas que deram
origem a todos os combustíveis utilizados para simulação do motor. Esses combustíveis
foram etanol anidro, óleo diesel e biodiesel. O software utilizado neste trabalho foi o
BOOST da empresa AVL. Para caracterizar os combustíveis, os seguintes dados são
requeridos pelo software:
Massa molar [kg/mol];
Poder Calorífico Inferior (PCI) [J/kg];
Razão ar/combustível estequiométrica;
Razão de massa de Carbono/massa total;
Razão de massa de Oxigênio/massa total;
Razão de massa de Nitrogênio/massa total;
Calor específico molar a pressão constante em função da temperatura
(𝑐�̅�) [kJ/kmol];
Entalpia específica molar em função da temperatura (ℎ̅) [kJ/kmol];
Entropia específica molar em função da temperatura (�̅�) [kJ/kmol];
Os valores de 𝑐�̅� , ℎ̅ e �̅� são calculados a partir de funções polinomiais que são
apresentadas nas equações 2, 3 e 4 [33].
Figura 6: Foto do motor da MAN instalado no banco de provas do LMT [32]
18
𝑐�̅� = 𝑅(𝑎1 + 𝑎2 ∗ 𝑇 + 𝑎3 ∗ 𝑇2 + 𝑎4 ∗ 𝑇3 + 𝑎5 ∗ 𝑇4) (2)
ℎ̅ = 𝑅(𝑎1 ∗ 𝑇 +𝑎2
2∗ 𝑇2 +
𝑎3
3∗ 𝑇3 +
𝑎4
4∗ 𝑇4 +
𝑎5
5∗ 𝑇5 + 𝑎6) (3)
�̅� = 𝑅(𝑎1 ∗ ln 𝑇 + 𝑎2 ∗ 𝑇 +𝑎3
2∗ 𝑇2 +
𝑎4
3∗ 𝑇3 +
𝑎5
4∗ 𝑇4 + 𝑎7) (4)
onde 𝑇 é a temperatura em Kelvin, 𝑅 é a constante universal dos gases em base molar, e
𝑎1 a 𝑎7 são constantes diferentes para cada combustível. Há dois valores para cada
conjunto de constantes de 𝑎1 a 𝑎7, sendo um substituído nas equações acima quando
T<1000 K (baixas temperaturas), e o outro, quando 𝑇 ≥ 1000 K (altas temperaturas). No
ponto de T=1000K os valores das equações acima, com as constantes de baixa e alta
temperaturas, são iguais. As equações não são válidas para qualquer temperatura, elas
têm limites superiores e inferiores. No BOOST, esses valores limites são diferentes para
cada substância. Os limites inferiores variam entre 200 K e 300 K, enquanto os limites
superiores variam entre 3500 K e 5000 K.
Para definir um combustível basta alterar os dados no arquivo “therm.dat” do
software BOOST. Os parâmetros massa molar [kg/mol], razão ar/combustível
estequiométrica, razão de Carbono/massa total, razão de Oxigênio/massa total e razão de
Nitrogênio/massa total, serão calculados automaticamente pelo programa quando a
molécula representativa da substância for definida, ou seja, quando for definido o número
de átomos de Carbono, Hidrogênio, Nitrogênio e Oxigênio que compõem uma molécula
representativa de uma dada substância. Vale ressaltar que o programa não aceita números
decimais para representar o número de átomos da molécula representativa.
O valor da constante 𝑎6 é definida a fim de deslocar os pontos de ℎ̅ na vertical.
Como essa temperatura é a temperatura base utilizada para calcular a entalpia de
formação, então o valor de 𝑎6 esta diretamente ligado ao valor de PCI.
A razão ar/combustível estequiométrica – (𝐴 𝐶⁄ )𝑒𝑠𝑡 – é uma razão entre a
quantidade mássica de ar mínima teórica necessária para que haja combustão completa
do combustível, ou seja, para que os produtos da combustão sejam apenas os gases inertes
presentes no ar atmosférico, mais 𝐻2𝑂 e 𝐶𝑂2 . De acordo com HEYWOOD [17], a
composição do ar atmosférico é aproximadamente a apresentada na tabela 3. Contudo,
como pode-se observar na quarta coluna, os gases inertes Ar e CO2 são considerados
como se fossem N2.
19
Tabela 3: composição do ar atmosférico [17].
Gás ppm de
volume
Massa
molecular
Fração
molar
Razão
molar
𝑶𝟐 209.500 31,998 0,2095 1
𝑵𝟐 780.900 28,012 0,7905 3,773
𝑨𝒓 9.300 39,348
𝑪𝑶𝟐 300 44,009
𝑨𝒓 𝒂𝒕𝒎𝒐𝒔𝒇é𝒓𝒊𝒄𝒐 1.000.000 28,962 1,0000 4,773
Para calcular (𝐴 𝐶⁄ )𝑒𝑠𝑡 dos combustíveis utilizados neste trabalho foi utilizado
o procedimento a seguir:
A. A fórmula geral para reação química de combustão completa pode ser
representada pela equação 5, sendo que Κ𝑂2, Κ𝐶𝑂2
, Κ𝐻2𝑂,𝑃 𝑒 Κ𝑂2 deverão ser
calculados pelas equações 6 a 9 e 𝑙, 𝑚, 𝑛 𝑒 Κ𝐻2𝑂,𝑅 são dados dos
combustíveis.
𝐶𝑙𝐻𝑚𝑂𝑛 + Κ𝑂2(𝑂2 + 3,773 ∗ 𝑁2) + Κ𝐻2𝑂,𝑅 ∗ 𝐻2𝑂
→ Κ𝐶𝑂2∗ 𝐶𝑂2 + Κ𝐻2𝑂,𝑃 ∗ 𝐻2𝑂 + Κ𝑁2
𝑁2 (5)
Κ𝐶𝑂2= 𝑙 (6)
Κ𝐻2𝑂,𝑃 =𝑚
2+ Κ𝐻2𝑂,𝑅; (7)
Κ𝑂2= 𝑙 +
𝑚
4−
𝑛
2; (8)
Κ𝑁2= Κ𝑂2
∗ 3,773; (9)
B. Sendo �̅�𝐶𝑜𝑚𝑏. (peso molecular do combustível) calculado na equação 10, e
�̅�𝐴𝑟 = 28,962 𝑘𝑔/𝑘𝑚𝑜𝑙 o peso molecular do ar, pode-se calcular (𝐴 𝐶⁄ )𝑒𝑠𝑡
pela equação 11.
�̅�𝐶𝑜𝑚𝑏. = 𝑙 ∗ �̅�𝐶 + 𝑚 ∗ �̅�𝐻 + 𝑛 ∗ �̅�𝑂 (10)
(𝐴 𝐶⁄ )𝑒𝑠𝑡 =Κ𝑂2
∗ (1 + 3,773) ∗ �̅�𝐴𝑟
1 ∗ �̅�𝐶𝑜𝑚𝑏. [𝑘𝑔 𝑘𝑔⁄ ] (11)
20
sendo �̅�𝐶 , �̅�𝐻 e �̅�𝑂 respectivamente o peso molecular dos átomos de Carbono,
Hidrogênio e Oxigênio. Seus valores são apresentados na tabela 4 [34, 35 e 36].
Tabela 4: Peso molecular de alguns átomos [34, 35 e 36].
Massa molecular kg/kmol
Carbono (�̅�𝑪) 12,0107
Oxigênio (�̅�𝑯) 1,00794
Hidrogênio (�̅�𝑶) 15,9994
4.2.1. Diesel
Buscou-se aproximar as propriedades do óleo diesel para o comercializado no
Brasil. Sendo assim, as propriedades que são apresentadas na tabela 5 foram calculadas a
partir dos dados de óleo diesel (sem biodiesel) encontrados no trabalho de FILHO [37].
Tabela 5: Propriedades do óleo diesel [37].
Propriedades Óleo diesel
Razão de Carbono/massa total 0,84615
Razão de Oxigênio/massa total 0
Razão de Nitrogênio/massa total 0
Poder Calorífico inferior [MJ/kg] 42,450
O BOOST já oferece um modelo de óleo diesel, e esses dados são apresentados
na tabela 6.
21
Tabela 6: Parâmetros do óleo diesel calculados pelo BOOST
Parâmetros Óleo diesel
Razão de Carbono/massa total 0,83905565
Razão de Oxigênio/massa total 0
Razão de Nitrogênio/massa total 0
Massa molar [kg/mol] 0,10020557
Poder Calorífico Inferior (PCI) [MJ/kg] 42,8390557
Razão ar/combustível estequiométrica – (𝑨 𝑪⁄ )𝒆𝒔𝒕 15,16942
Como pode-se observar comparando os dados da tabela 6 e 5, há uma pequena
diferença entre as razão de Carbono/massa total, porém, como essa diferença é pequena
decidiu-se manter a molécula de óleo diesel apresentada pelo BOOST. A única exceção
foi o valor de PCI que será ajustado a fim de igualar o apresentado por FILHO.
O BOOST também apresenta os dados termodinâmicos do óleo diesel. Portanto,
as constantes 𝑐1 − 𝑐5 e 𝑐7 de 𝑐�̅� , ℎ̅ e �̅� serão aproveitadas, porém o valor de 𝑎6 para
baixas temperaturas foi modificado a fim de ajustar o PCI. O parâmetro 𝑎6 para altas
temperaturas também foi modificado, entretanto, apenas para que, em 𝑇 = 1000 𝐾, as
curvas da entalpia para altas e baixas temperaturas sejam equivalentes.
Apesar de 𝑎6 ser alterado no arquivo “therm.dat” para ajustar o PCI do
combustível, o programa sempre apresenta o valor de 𝑎6 ajustado para ℎ̅ = 0 em 25ºC.
Os valores das constantes apresentadas pelo programa ao usuário são apresentadas na
tabela 7.
22
Tabela 7: Constantes para determinar as características termodinâmicas do óleo diesel
Constantes Baixas
temperaturas
Altas
temperaturas
𝒂𝟏 -2,76665E+00 1,89304E+01
𝒂𝟐 9,82799E-02 4,87857E-02
𝒂𝟑 -3,49532E-05 -1,68215E-05
𝒂𝟒 -2,29927E-08 2,69602E-09
𝒂𝟓 1,59497E-11 -1,65882E-13
𝒂𝟔 -3,19659E+03 -9,38951E+12
𝒂𝟕 3,51970E+01 -7,88116E+01
4.2.2. Etanol
Assim como no caso do óleo diesel, foi encontrado no trabalho de MELO [5] as
propriedades de um etanol hidratado, que são apresentadas na tabela 8.
Tabela 8: Propriedades do etanol hidratado [5].
Parâmetros Etanol
Massa específica (𝒎𝒆𝒕𝒂𝒏𝒐𝒍𝑯) [kg/m3] 808,7
Poder Calorífico Inferior (𝑷𝑪𝑰𝒆𝒕𝒂𝒏𝒐𝒍𝑯) [MJ/kg] 24,76
Razão ar/combustível estequiométrica – (𝑨 𝑪⁄ )𝒆𝒔𝒕 8,8
Razão de Carbono/massa total (𝑹𝑪) 0,507
Razão de Oxigênio/massa total (𝑹𝑶) 0,363
Razão de Nitrogênio/massa total (𝑹𝑵) 0
Etanol (E) [% v/v] 95,7
Água (A) [% v/v] 4,3
Neste trabalho deseja-se realizar misturas utilizando etanol anidro. De acordo
com a resolução Nº 7 de 2011 da ANP [38], o etanol anidro pode conter no máximo 0,4%
v/v de água. Será considerado neste trabalho etanol puro, ou seja, 0% v/v de água. Sendo
23
assim, para encontrar os valores necessários para alimentar o BOOST, foram realizados
os cálculos que serão apresentados a seguir.
Considerando que a massa específica da água a 1atm e 25ºC (𝑚𝐻2𝑂) é 997 kg/m3
[39] e que a sua massa molecular (�̅�𝐻2𝑂) é 18,015 kg/kmol [39] então, pode-se calcular
a massa específica de etanol puro (𝑚𝑒𝑡𝑎𝑛𝑜𝑙𝑃) pela equação 12.
𝑚𝑒𝑡𝑎𝑛𝑜𝑙𝑃 =𝑚𝑒𝑡𝑎𝑛𝑜𝑙𝐻 − (A 100⁄ ) ∙ 𝑚𝐻2𝑂
(E 100⁄ ) (12)
Para encontrar a razão de Carbono e Hidrogênio em relação à massa total do
etanol puro, primeiro devemos calcular a massa específica de cada átomo na mistura de
etanol hidratado apresentado nas equações de 13 a 15.
𝑚𝐶,𝑒𝑡𝑎𝑛𝑜𝑙𝐻 = 𝑚𝑒𝑡𝑎𝑛𝑜𝑙𝐻 ∙ 𝑅𝐶 (13)
𝑚𝑂,𝑒𝑡𝑎𝑛𝑜𝑙𝐻 = 𝑚𝑒𝑡𝑎𝑛𝑜𝑙𝐻 ∙ 𝑅𝑂 (14)
𝑚𝐻,𝑒𝑡𝑎𝑛𝑜𝑙𝐻 = 𝑚𝑒𝑡𝑎𝑛𝑜𝑙𝐻 ∙ [1 − (𝑅𝐶 + 𝑅𝑂)] (15)
Sendo 𝑚𝐶,𝑒𝑡𝑎𝑛𝑜𝑙𝐻 a massa específica de Carbono no etanol hidratado,
𝑚𝑂,𝑒𝑡𝑎𝑛𝑜𝑙𝐻 a de Oxigênio e 𝑚𝐻,𝑒𝑡𝑎𝑛𝑜𝑙𝐻 a de Hidrogênio.
As massas específicas do carbono e do hidrogênio da água são calculadas pelas
equações 16 e 17.
𝑚𝑂, 𝐻2𝑂 =𝑚𝐻2𝑂
�̅�𝐻2𝑂∙ 1 ∙ �̅�𝑂 ∙ A (16)
𝑚𝐻, 𝐻2𝑂 =𝑚𝐻2𝑂
�̅�𝐻2𝑂∙ 2 ∙ �̅�𝐻 ∙ A (17)
Podemos então calcular a razão mássica de cada componente para o etanol puro
(etanolP) utilizando as equações 18 a 20.
𝑅𝐶,𝑒𝑡𝑎𝑛𝑜𝑙𝑃 =𝑚𝐶,𝑒𝑡𝑎𝑛𝑜𝑙𝐻
𝑚𝑒𝑡𝑎𝑛𝑜𝑙𝑃 ∙ (E 100⁄ ) (18)
𝑅𝑂,𝑒𝑡𝑎𝑛𝑜𝑙𝑃 =𝑚𝑂,𝑒𝑡𝑎𝑛𝑜𝑙𝐻 − 𝑚𝑂, 𝐻2𝑂
𝑚𝑒𝑡𝑎𝑛𝑜𝑙𝑃 ∙ (E 100⁄ ) (19)
𝑅𝐻,𝑒𝑡𝑎𝑛𝑜𝑙𝑃 =𝑚𝐻,𝑒𝑡𝑎𝑛𝑜𝑙𝐻 − 𝑚𝐻, 𝐻2𝑂
𝑚𝑒𝑡𝑎𝑛𝑜𝑙𝑃 ∙ (E 100⁄ ) (20)
24
Calculando as equações acima temos, na tabela 9, as razões mássicas de carbono,
oxigênio e Hidrogênio em relação à massa total de etanol com 0 % v/v de H2O.
Tabela 9: Razões mássicas do etanol puro (0 % v/v de H2O)
Elemento químico Razão mássica
Razão de Carbono/massa total (𝑹𝑪,𝒆𝒕𝒂𝒏𝒐𝒍𝑷) 0,535
Razão de Oxigênio/massa total (𝑹𝑶,𝒆𝒕𝒂𝒏𝒐𝒍𝑷) 0,334
Razão de Hidrogênio/massa total (𝑹𝑯,𝒆𝒕𝒂𝒏𝒐𝒍𝑷) 0,131
De acordo com VLASSOV [40] pode-se aproximar a molécula do etanol
brasileiro para C2H5OH. Nota-se então, que os valores das razões mássicas de cada átomo
calculados para a molécula C2H5OH, apresentados na tabela 10, são semelhantes aos
encontrados na tabela 9.
Tabela 10: Razões mássicas do etanol da molécula C2H5OH
Elemento químico Razão mássica
Carbono 0,521
Oxigênio 0,131
Hidrogênio 0,347
O Poder Calorífico Inferior (PCI) de um combustível é a energia liberada,
geralmente à pressão constante [17], quando se tem combustão completa de uma
substância e as moléculas de H2O estão no estado de vapor. Sendo assim, para calcular o
PCI do etanol puro (0% de água nos reagentes) a partir do PCI do etanol hidratado,
calculado por MELO [5], deve-se computar a diferença entre a entalpia de formação do
H2O no estado liquido (nos reagentes) e no estado de vapor (nos produtos). De acordo
com BORGNAKKE [39], as entalpias de formação da água liquida (∆ℎ̅𝑓,𝑙0 ) e da água em
estado de vapor (∆ℎ̅𝑓,𝑔0 ) possuem os valores apresentados na tabela 11.
25
Tabela 11: Entalpias de formação da água liquida (l) e em forma de vapor (g)
∆�̅�𝒇,𝒍𝟎 -241.826 kJ/kmol
∆�̅�𝒇,𝒈𝟎 -285.830 kJ/kmol
O PCI do etanol puro é calculado pelas equações 21 e 22, se baseando nos valores
apresentados na tabela 11.
∆ℎ𝑓,𝑙→𝑔 =∆ℎ̅𝑓,𝑔
0 − (∆ℎ̅𝑓,𝑔0 )
�̅�𝐻2𝑂 [
𝑘𝐽
𝑘𝑔] (21)
𝑃𝐶𝐼𝑒𝑡𝑎𝑛𝑜𝑙𝑃 = [𝑃𝐶𝐼𝑒𝑡𝑎𝑛𝑜𝑙𝐻 −𝐴
𝑚𝑒𝑡𝑎𝑛𝑜𝑙𝐻∙ 𝑚𝐻2𝑂 ∙ ∆ℎ𝑓,𝑙→𝑔] ∗
𝑚𝑒𝑡𝑎𝑛𝑜𝑙𝐻
𝑚𝑒𝑡𝑎𝑛𝑜𝑙𝑃 (22)
Calculando a equação acima e (𝐴 𝐶⁄ )𝑒𝑠𝑡 (descrito em 4.2), utilizando C2H5OH
como molécula representativa do etanol, foi montada a tabela 12 com dados do etanol
puro.
Tabela 12: Propriedades do etanol puro
Parâmetros Etanol
Poder Calorífico Inferior (PCI) [MJ/kg]; 25,1526
Massa molar [kg/mol]; 0,0461
Razão ar/combustível estequiométrica – (𝑨 𝑪⁄ )𝒆𝒔𝒕; 9,002
Assim como no caso do óleo diesel, o etanol já era um combustível existente no
BOOST. Sendo assim, as variáveis 𝑎1 − 𝑎5 e 𝑎7 também foram aproveitadas, porém,
como a molécula padrão foi alterada essas constantes foram ajustadas a fim de manter os
valores de 𝑐�̅� , ℎ̅ e �̅� na base mássica iguais aos previamente existentes no programa.
Sendo assim, as constantes foram multiplicadas por: �̅�𝑒𝑡𝑎𝑛𝑜𝑙,𝑛𝑜𝑣𝑜 �̅�𝑒𝑡𝑎𝑛𝑜𝑙,𝑎𝑛𝑡𝑖𝑔𝑜⁄ . Já as
constantes 𝑎6 foram ajustadas a fim de se obter o valor de PCI desejado e manter a
continuidade da curva de ℎ̅, analogamente ao que foi realizado no caso do óleo diesel. Os
valores apresentados pelo programa se encontram na tabela 13.
26
Tabela 13: Constantes para determinar as características termodinâmicas do etanol
puro.
Constantes Baixas
temperaturas
Altas
temperaturas
𝒂𝟏 0,48586957E+01 0,65624365E+01
𝒂𝟐 -0,37401726E-02 0,15204222E-01
𝒂𝟑 0,69555378E-04 -0,53896795E-05
𝒂𝟒 -0,88654796E-07 0,86225011E-09
𝒂𝟓 0,35168835E-10 -0,51289787E-13
𝒂𝟔 -1,7383035E+04 -3,26779231E+03
𝒂𝟕 0,48018545E+01 -0,94730202E+01
4.2.3. Biodiesel
Buscando se aproximar das propriedades do biodiesel produzido no Brasil, foi
feita uma pesquisa bibliográfica dos dados de biodiesel a base de soja, que é a matéria
prima principal na produção de biodiesel no país.
LIN et al. [41] estudaram alguns biodieseis, sendo que entre eles se encontrava
uma amostra de biodiesel comercial produzida a partir de óleo de soja e álcool metílico
via reação de transesterificação (análogo ao processo de produção predominante no
Brasil). As propriedades deste combustível são apresentadas na tabela 14.
Tabela 14: Propriedades do biodiesel encontrado no trabalho de LIN et al. [41].
Parâmetros Biodiesel
Poder Calorífico Inferior (PCI) [MJ/kg]; 39,714
Razão de Carbono/massa total; 0,7696
Razão de Oxigênio/massa total; 0,0941
Razão de Nitrogênio/massa total; 0,0000
27
Para definir a molécula padrão do biodiesel de soja, algumas considerações
foram feitas. LIU et al. [42] realizaram a análise de uma amostra de biodiesel a base de
óleo de soja e álcool metílico, encontrando a razão mássica de diferentes moléculas
presentes no combustível. Os resultados são apresentados na tabela 15.
Tabela 15: Análise da composição do biodiesel.
Molécula Nº de
Carbonos (C)
Percentual
mássico
Percentual
mássico por C
C14:0 14 0
C16:0 16
10,81 10,92
C16:1 0,11
C18:0
18
4,54
87,43 C18:1 24,96
C18:2 50,66
C18:3 7,27
C20:0 20
0,37
0,69 C20:1 0,32
C22:0 22 0,42 0,42
C24:0 24 0,12 0,12
Total 99,58 99,58
Como pode ser observado, a maioria das moléculas que compõem o biodiesel à
base de soja são compostos por 18 carbonos. Sendo assim, a molécula padrão escolhida
terá 18 carbonos. Buscando aproximar-se ao máximo das razões mássicas de carbono,
hidrogênio e oxigênio apresentadas na tabela 14, foi definida a molécula padrão do
biodiesel a base de soja como C18H33O2. Logo, as novas razões mássicas de cada átomo
presente no combustível são apresentadas na tabela 16, enquanto o PCI será mantido o
apresentado na tabela 14. Na tabela 16, também são apresentado os valores de (𝐴 𝐶⁄ )𝑒𝑠𝑡
e Massa Molar, calculados como descrito no item 4.2.
28
Tabela 16: Propriedades do biodiesel, parte 2 de 2
Parâmetros Biodiesel
Poder Calorífico Inferior (PCI) [MJ/kg]; 39,714
Razão de Carbono/massa total 0,7681
Razão de Hidrogênio/massa total 0,1182
Razão de Oxigênio /massa total; 0,1137
Razão ar/combustível estequiométrica – (𝑨 𝑪⁄ )𝒆𝒔𝒕 12,40
Massa molar [kg/mol] 0,2815
O biodiesel não foi encontrado na base de dados do BOOST. Sendo assim, as
propriedades termodinâmicas do biodiesel serão aproximadas para a do diesel. Como 𝑐�̅�,
ℎ̅ e �̅� se encontram na base molar, para que seus valores na base mássica sejam
equivalentes ao do óleo Diesel, as constantes serão multiplicadas pela razão da massa
molar do biodiesel com o diesel (�̅�𝑏𝑖𝑜𝑑𝑖𝑒𝑠𝑒𝑙 �̅�𝑑𝑖𝑒𝑠𝑒𝑙⁄ ). Em seguida, a constante 𝑎6 foi
ajustada a fim de ajustar o PCI, e 𝑎6 para altas temperaturas para que, em 𝑇 = 1000 𝐾,
as curvas da entalpia para altas e baixas temperaturas sejam equivalentes. Os valores das
constantes 𝑎1 − 𝑎7 para o biodiesel se encontram na tabela 17, lembrando que o programa
BOOST apresenta o valor de 𝑎6 ajustado para ℎ̅ = 0 em 25ºC
Tabela 17: Constantes para determinar as características termodinâmicas do Biodiesel
Constantes Baixas
temperaturas
Altas
temperaturas
𝒂𝟏 -7,77111E+00 5,31727E+01
𝒂𝟐 2,76054E-01 1,37032E-01
𝒂𝟑 -9,81783E-05 -4,72491E-05
𝒂𝟒 -6,45831E-08 7,57272E-09
𝒂𝟓 4,45194E-11 -4,65938E-13
𝒂𝟔 -2,81752E+05 -8,97877E+09
𝒂𝟕 9,88631E+01 -2,21370E+02
29
4.2.4. Misturas de combustíveis
Esse trabalho analisará misturas de combustíveis envolvendo óleo diesel,
biodiesel e etanol. Como pode ser observado nos trabalhos de AL-HASSAN et al [9] e
ROCHA [7], misturas de diesel, biodiesel e etanol com até 5% v/v de etanol são possíveis
de se manter estáveis.
A Medida Provisória 647/2014 [43], que está sobre análise do poder legislativo
brasileiro, prevê aumento no teor de biodiesel no óleo diesel para 6% a partir de Julho e
7% a partir de Novembro de 2014.
Nos EUA, quase todas as montadoras permitem o uso de óleo diesel com até
20% v/v de biodiesel sem que a garantia de fábrica do veículo seja comprometida [44].
Levando em conta os fatos descritos acima, as misturas usadas neste trabalho
foram escolhidas e são apresentadas na tabela 18. A sigla E representa etanol e o número
a sua direita o %v/v de etanol na mistura. Analogamente, a sigla B representa biodiesel e
o número a sua direita o %v/v de biodiesel na mistura. O valor restante necessário para
completar 100% é o %v/v de óleo diesel na mistura.
Tabela 18: Misturas de combustíveis escolhidas para simulação
Etanol
Biodiesel
1% v/v 2% v/v 5% v/v
7% v/v E1B7 E2B7 E5B7
10% v/v E1B10 E2B10 E5B10
20% v/v E1B20 E2B20 E5B20
4.3. Modelo de combustão
Utilizando a 1ª Lei da Termodinâmica e os limites da câmara de combustão como
condições de contorno, pode-se calcular a variação de energia interna no interior do
cilindro. Essa variação de energia envolve o trabalho no pistão, o calor liberado pelo
combustível, a transferência de calor pelas paredes, vazão de entalpia devido ao “blow-
by” (passagem de massa entre os anéis do pistão e a parede da câmara de combustão) e a
30
entalpia presente na massa que entra e sai do cilindro pelas válvulas de admissão e
exaustão [39].
De acordo com HEYWOOD [17], a equação 23, baseada na 1ª Lei da
Termodinâmica, representa a taxa de variação de energia interna do cilindro em função
do ângulo do virabrequim.
𝑑(𝑚 ∙ 𝑢)
𝑑𝜃= −𝑃𝑐 ∙
𝑑𝑉
𝑑𝜃+
𝑑𝑄𝑐𝑜𝑚𝑏
𝑑𝜃− ∑
𝑑𝑄𝑤
𝑑𝜃− ℎ𝐵𝐵 ∙
𝑑𝑚𝐵𝐵
𝑑𝜃+ ∑
𝑑𝑚𝑒
𝑑𝜃∙ ℎ𝑒
− ∑𝑑𝑚𝑠
𝑑𝜃∙ ℎ𝑠 − 𝑞𝑒𝑣 ∙ 𝑓 ∙
𝑑𝑚𝑒𝑣
𝑑𝜃
(23)
Em que:
𝜃 é o ângulo do virabrequim;
𝑑(𝑚𝑐∙𝑢)
𝑑𝜃 é a taxa de variação de energia interna do cilindro em função de 𝜃;
𝑃𝑐 ∙𝑑𝑉
𝑑𝜃 é o trabalho realizado pelo pistão em função de 𝜃;
𝑑𝑄𝑐𝑜𝑚𝑏
𝑑𝜃 é a taxa de liberação de calor pelo combustível em função de 𝜃;
∑𝑑𝑄𝑤
𝑑𝜃 é a taxa de calor perdido pela parede em função de 𝜃;
ℎ𝐵𝐵 ∙𝑑𝑚𝐵𝐵
𝑑𝜃 é o fluxo de entalpia devido ao blow-by em função de 𝜃;
𝑑𝑚𝑒 é o elemento de massa entrando no cilindro;
ℎ𝑒 é a entalpia do elemento de massa entrando no cilindro;
𝑑𝑚𝑠 é o elemento de massa saíndo no cilindro;
ℎ𝑠 é a entalpia do elemento de massa saindo no cilindro;
𝑞𝑒𝑣 é o calor de vaporização do combustível;
𝑓 é o fator de vaporização do combustível;
𝑚𝑒𝑣 é a massa de combustível evaporada.
Uma representação simplificada do balanço de energia no cilindro do motor é
apresentada pelo esquema na figura 7.
31
Figura 7: Esquema do balanço de energia no motor
Foi aplicado neste trabalho o modelo de combustão de Wiebe duas zonas para
simular a combustão. O modelo de combustão de Wiebe determina a taxa de liberação de
calor da combustão, porém, antes da explicação do modelo de combustão de Wiebe, as
equações das duas zonas serão apresentadas.
As duas zonas correspondem à massa queimada e à massa não queimada no
interior do cilindro. Considera-se que a pressão das duas zonas são iguais porém a
temperatura é diferente. Baseando-se na 1ª Lei da Termodinâmica, as equações 24 e 25
representam a taxa de variação de energia interna do cilindro, em função do ângulo do
virabrequim, da massa queimada e da massa não queimada.
𝑑(𝑚𝑏 ∙ 𝑢𝑏)
𝑑𝜃= −𝑃𝑐 ∙
𝑑𝑉𝑏
𝑑𝜃+
𝑑𝑄𝑐𝑜𝑚
𝑑𝜃− ∑
𝑑𝑄𝑤,𝑏
𝑑𝜃+
𝑑𝑚𝑏−𝑢
𝑑𝜃∙ ℎ𝑏−𝑢
−𝑑𝑚𝐵𝐵,𝑏
𝑑𝜃∙ ℎ𝐵𝐵,𝑏
(24)
𝑑(𝑚𝑢 ∙ 𝑢𝑢)
𝑑𝜃= −𝑃𝑐 ∙
𝑑𝑉𝑢
𝑑𝜃− ∑
𝑑𝑄𝑤,𝑢
𝑑𝜃+
𝑑𝑚𝑏−𝑢𝑛𝑏
𝑑𝜃∙ ℎ𝑏−𝑢𝑛𝑏
−𝑑𝑚𝐵𝐵,𝑢
𝑑𝜃∙ ℎ𝐵𝐵,𝑢
(25)
32
onde:
O sufixo 𝑏 indica massa queimada e o sufixo 𝑢 indica massa não queimada.
𝑑𝑚𝑏−𝑢𝑛𝑏
𝑑𝜃∙ ℎ𝑏−𝑢𝑛𝑏 é o termo responsável pelo cálculo da transferência de entalpia
entre a massa queimada e não queimada devido à combustão de um determinado
elemento de massa. A transferência de calor entre as massas queimadas e não
queimadas é desprezada neste modelo.
4.3.1. Modelo de combustão de Wiebe
A taxa de queima da mistura ar-combustível, o início da combustão e sua
duração são influenciados por muitos parâmetros, tais como a geometria do motor, a
turbulência no interior do cilindro, temperaturas, características do combustível, massa
de ar admitida, etc. Há inúmeras variáveis que influenciam à combustão.
O modelo de combustão de Wiebe [17], utilizado neste trabalho, utiliza a
equação 26 para prever a fração de massa queimada de combustível em função do ângulo
do virabrequim. Este modelo considera inicialmente o início e a duração da combustão.
Em seguida, utilizando os parâmetros 𝑎 e 𝑚 , a equação busca considerar as outras
variáveis que influenciam a queima da massa de combustível.
𝑥(𝜃) = 1 − exp [−𝑎 ∗ (
𝜃 − 𝜃𝑖
∆𝜃)
𝑚+1
] (26)
Na equação acima, 𝜃 é o ângulo do virabrequim, 𝑥(𝜃) é a fração de massa
queimada até o ângulo 𝜃, 𝜃𝑖 é o ângulo do início da liberação de energia (ou início da
combustão), ∆𝜃 é a duração da combustão e 𝑎 e 𝑚 são parâmetros a serem ajustados
junto da curva experimental de liberação de calor.
Pode ser observado pela equação 26 que o parâmetro 𝑎 está diretamente
relacionado à massa total de combustível queimada no intervalo ∆𝜃. No momento que 𝜃
é o ângulo do fim da combustão, (𝜃 − 𝜃𝑖) = ∆𝜃, então a equação 26 pode ser simplificada
pela equação 27, no final da combustão.
𝑥(𝜃) = 1 − exp [−𝑎] (27)
33
Sendo assim, o parâmetro 𝑎 definirá a massa total de combustível queimada
durante o intervalo considerado por Wiebe como a duração da combustão ∆𝜃.
Na figura 8 pode-se observar uma curva de Wiebe generalizada.
Figura 8: Modelo geral da fração mássica de combustível queimada do modelo de
Wiebe
Foi feita uma pesquisa bibliográfica buscando encontrar valores dos parâmetros
𝑎 e 𝑚 para as misturas de combustíveis utilizadas. Apesar de consultar diversos bancos
de dados, não houve êxito em encontrar valores desses parâmetros. No entanto, no
trabalho de GOGOI et al. [45] foram encontrada curvas da taxas de liberação de calor,
figura 9, dos seguintes combustíveis: diesel (E0B0), biodiesel de Pinhão-manso (Jatropha
curcas) (E0B100) e uma mistura de 80% de biodiesel com 20% de etanol (E20B80).
0
0,2
0,4
0,6
0,8
1
0 5 10 15 20 25 30 35 40 45
∆𝜃
𝜃𝑖 𝜃𝑓
𝑥(𝜃
)
Ângulo
34
Figura 9: Taxa liquida de liberação de calor de 3 combustíveis em um dado motor
Utilizando o software FindGraph, as curvas de E0B0, E0B100 e E20B80, com
taxas de liberação de calor positivas, foram convertidas para dados numéricos. Em
seguida esses dados foram integrados em função do ângulo do virabrequim, a fim de
encontrar a curva da fração de calor liberada em função do ângulo do virabrequim.
A curva de Wiebe foi ajustada para a fração de massa queimada entre 1% e 99%
das curvas da figura 9. Foi estipulado 𝑎 = 4,6, o que representa aproximadamente 99%
de massa queimada pela equação de Wiebe. O parâmetro 𝑚 foi ajustado manualmente,
equiparando um ponto arbitrário da curva experimental com a curva teórica. O ponto
escolhido variou para cada curva a fim de encontrar a melhor aproximação da curva de
Wiebe com a curva experimental. As figuras 10, 11 e 12 apresentam respectivamente as
curvas de Wiebe ajustadas (linha tracejada) e as curvas teóricas (linha contínua) de E0B0,
E0B100 e E20B80.
35
Figura 10: Fração de massa queimada experimental e ajuste da curva de Wiebe para
óleo diesel (E0B0).
Figura 11: Fração de massa queimada experimental e ajuste da curva de Wiebe para
biodiesel (E0B100).
36
Figura 12: Fração de massa queimada experimental e ajuste da curva de Wiebe para
E20B80.
Os valores de 𝑚 para cada combustível são encontrados na tabela 19:
Tabela 19: parâmetro m ajustado para E0B0, E0B100 e E20B80
Combustível 𝒎
E0B0 (𝒎𝒅) 0,7894
E0B100 (𝒎𝒃) 0,7747
E20B80 (𝒎𝒃𝒆) 0,7792
Como pode-se observar, não há grande diferença entre os valores do parâmetro
𝑚 em função do combustível. Sendo assim, foi suposto que seria uma boa aproximação
considerar que cada combustível influencia linearmente o valor de 𝑚 para cada mistura.
Primeiro foi identificado o valor de 𝑚 para o etanol puro (𝑚𝑒). Em seguida, foi calculado
o valor de 𝑚 para cada mistura que será utilizada neste trabalho. Tais valores são
apresentados na tabela 20, sendo o cálculo realizado utilizando a equação 28 para cada
mistura.
37
𝑚 =
𝐷 ∙ 𝑚𝑑 + 𝐵 ∙ 𝑚𝑏 + 𝐸 ∙ 𝑚𝑒
100 (28)
onde D, B e E são respectivamente os percentuais de óleo diesel, biodiesel e etanol na
mistura.
Tabela 20: Valores calculados do parâmetro m do modelo de Wiebe
Combustível 𝒎
E1B7 0,7885
E1B10 0,7880
E1B20 0,7866
E2B6 0,7886
E2B10 0,7881
E2B20 0,7866
E5B6 0,7795
E5B10 0,7883
E5B20 0,7869
4.4. Modelo de transferência de calor
A transmissão de calor entre o gás no interior do cilindro e as paredes do cilindro
é calculada pela equação 29. Nesta equação, 𝑄𝑖 é o fluxo de calor transmitido para a
parede, 𝐴𝑖 é a área da superfície, ℎ é o coeficiente de transferência de calor, 𝑇𝑐 é a
temperatura média do gás no cilindro e 𝑇𝑖 é a temperatura da parede. O ℎ é de difícil
determinação teórica. Sendo assim, este é estimado a partir de modelos.
𝑄𝑖 = ℎ ∙ 𝐴𝑖 ∙ (𝑇𝑐 − 𝑇𝑖) (29)
Para a simulação do motor, foi utilizado o modelo de transmissão de calor criado
pela empesa AVL [14], que será representado pelo símbolo ℎ𝐴𝑉𝐿. Este modelo se baseia
no modelo de Woschni, publicado em 1990 (ℎ𝑤,1990 ) [46], que por sua vez é um
aperfeiçoamento do modelo publicado por Woschni, em 1967 (ℎ𝑤,1967) [19]. ℎ𝑤,1967 é
apresentado na equação 30.
38
ℎ𝑤,1967 = 3,26 ∙ 𝐷−0,2 ∙ 𝑃𝑐
0,8 ∙ 𝑇𝑐−0,55 [𝐶1 ∙ 𝑐𝑚 + 𝐶2(𝑃𝑐 − 𝑃0)
𝑉𝑑 ∙ 𝑇1
𝑃1 ∙ 𝑉1]
0,8
(30)
onde:
D é o diâmetro do cilindro;
𝑃𝑐 é a pressão média dos gases no interior do cilindro;
𝑇𝑐 é a temperatura média dos gases no interior do cilindro;
𝐶1 = {6,18 + 0,417 ∙ 𝑐𝑢 𝑐𝑚⁄ , 𝑝𝑒𝑟í𝑜𝑑𝑜 𝑑𝑒 𝑒𝑛𝑡𝑟𝑎𝑑𝑎 𝑒 𝑠𝑎í𝑑𝑎 𝑑𝑒 𝑓𝑙𝑢𝑖𝑑𝑜
2,28 + 0,308 ∙ 𝑐𝑢 𝑐𝑚⁄ , durante os outros períodos do motor
𝑐𝑢 = 𝐷 ∙ 𝜔𝑝 2⁄
𝜔𝑝 é a velocidade angular da pá usada para medir o turbilhonamento do gás;
𝑐𝑚 é a velocidade média do pistão;
𝐶2 = {0,00324, 𝑑𝑢𝑟𝑎𝑛𝑡𝑒 𝑎 𝑒𝑥𝑝𝑎𝑛𝑠ã𝑜 𝑒 𝑎 𝑐𝑜𝑚𝑏𝑢𝑠𝑡ã𝑜; 0, 𝑑𝑢𝑟𝑎𝑛𝑡𝑒 𝑎 𝑐𝑜𝑚𝑝𝑟𝑒𝑠𝑠ã𝑜, 𝑒 𝑎 𝑒𝑛𝑡𝑟𝑎𝑑𝑎 𝑒 𝑎 𝑠𝑎í𝑑𝑎 𝑑𝑒 𝑓𝑙𝑢𝑖𝑑𝑜;
𝑃0 é a pressão no interior da câmara de combustão sem ocorrência de combustão
(“motored”);
𝑉𝑑 é o volume deslocado;
𝑇1, 𝑃1 e 𝑉1 são respectivamente a temperatura, a pressão e o volume no ângulo de
fechamento da válvula de admissão.
O modelo de Woschni publicado em 1990 (ℎ𝑤,1990) é apresentado pela equação
31.
ℎ𝑤,1990 = 3,26 ∙ 𝐷−0,2 ∙ 𝑃𝑐
0,8 ∙ 𝑇𝑐−0,55 {𝐶1 ∙ 𝑐𝑚 [1 + 2 (
𝑉𝑐
𝑉)
2
𝐼𝑀𝐸𝑃−0,2]}
0,8
(31)
onde:
𝑉𝑐 é o volume morto do cilindro;
𝑉 é o volume instantâneo no interior do cilindro;
𝐼𝑀𝐸𝑃 é a pressão média efetiva indicada.
Apesar do modelo de transferência de calor de Woschni ℎ𝑤,1990 apresentar
melhores resultados que o modelo que o antecedeu, ele é valido para uma faixa mais
restrita de funcionamento do motor [47]. Por esse motivo, ele será utilizado sempre que
a equação 32 for verdadeira, caso contrário deve-se utilizar o modelo de Woschni
publicado em 1967 (ℎ𝑤,1967).
39
𝐶1 ∙ 𝑐𝑚 [1 + 2 (
𝑉𝑐
𝑉)
2
𝐼𝑀𝐸𝑃−0,2] > 𝐶2(𝑃𝑐 − 𝑃0)𝑉𝑑 ∙ 𝑇1
𝑃1 ∙ 𝑉1 (32)
O modelo de coeficiente de calor proposto pela empresa AVL [47] se baseia no
modelo de Woschni. Tal modelo é apresentado pela equação 33.
ℎ𝐴𝑉𝐿 = {
𝑐𝑜𝑝𝑟𝑒𝑠𝑠ã𝑜 𝑒 𝑒𝑥𝑝𝑎𝑚𝑠ã𝑜: ℎ𝑤
𝑎𝑑𝑚𝑖𝑠𝑠ã𝑜 𝑒 𝑒𝑥𝑎𝑢𝑠𝑡ã𝑜: {ℎ𝑤, 𝑠𝑒 ℎ𝑤 > ℎ𝑤 𝑚𝑜𝑑 ℎ𝑤 𝑚𝑜𝑑, 𝑠𝑒 ℎ𝑤 𝑚𝑜𝑑 > ℎ𝑤
(33)
onde:
ℎ𝑤 = {ℎ𝑤,1990, 𝑠𝑒 𝑎 𝑒𝑞𝑢𝑎çã𝑜 32 𝑓𝑜𝑟 𝑣𝑒𝑟𝑑𝑎𝑑𝑒;
ℎ𝑤,1967, 𝑐𝑎𝑠𝑜 𝑐𝑜𝑛𝑡𝑟á𝑟𝑖𝑜;
ℎ𝑤 𝑚𝑜𝑑 = 3,26 ∙ 𝐷−0,2 ∙ 𝑃𝑐0,8 ∙ 𝑇𝑐
−0,55 (14,0 (𝑑𝑖𝑛
𝐷)
2|𝑉𝑖𝑛|)
0,8
o 𝐷 é o diâmetro do cilindro;
o 𝑃𝑐 e 𝑇𝑐 são respectivamente a pressão e a temperatura médias do
gás no interior do cilindro;
o 𝑑𝑖𝑛 é o diâmetro do duto conectado à válvula de admissão;
o 𝑉𝑖𝑛 é a velocidade do gás na válvula de entrada;
4.5. Programa para análise das curvas de pressão experimentais
O Laboratório de Maquinas Térmicas (LMT) da UFRJ, onde o motor em que foi
desenvolvido este trabalho se encontra, cedeu curvas de pressão experimentais. Foram
três tomadas de duzentos ciclos da pressão no interior do cilindro para cada uma das
seguintes potências: 25%, 50%, 75% e 100% da potência nominal do motor. Os
programas que se encontram no anexo I e II foram escritos no software MATLab a fim
de calcular os seguintes dados: início da combustão, fim da combustão, potência indicada
(potência produzida levando em conta perda de calor, porém não computando perdas por
atrito) e potência dissipada pelo atrito (calculada levando em conta a potência nominal do
motor e a potência indicada). O anexo I utiliza o método descrito no item 4.5.2.1 para
calcular o fim da combustão, enquanto o anexo II utiliza o método descrito no 4.5.2.2.
40
4.5.1. Início da combustão
De acordo com HEYWOOD [17], tanto a compressão quanto a expansão,
quando não estão somadas ao fenômeno da combustão, apresentam uma forte correlação
com a equação 34.
𝑃𝑉𝑘 = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑒 (34)
𝑘 =𝑐𝑝
𝑐𝑣
(35)
em que 𝑃 e 𝑉 são respectivamente, a pressão e o volume do gás no interior do cilindro; e
𝑘, apresentado na equação 35, é a razão entre os calores específicos do gás à pressão
constante (𝑐𝑝) e à volume constante (𝑐𝑣). A combustão ocorre quando a equação 34 não
é obedecida entre as etapas de compressão e expansão do cilindro. Para identificar o início
da combustão, basta encontrar o momento em que há uma mudança brusca do
comportamento dos gases no interior do cilindro em relação à equação 34. Já no caso do
fim da combustão, essa mudança é bastante sutil, o que torna mais difícil identifica-la.
Considerando a pressão, foram utilizados os dados experimentais cedidos pelo
LMT, que se encontram em função do ângulo do virabrequim. O volume instantâneo no
interior do cilindro foi calculado teoricamente em função do ângulo do virabrequim,
utilizando a equação 36, apresentada por HEYWOOD [17].
𝑉(𝜃) = {1 +
1
2[𝑅 + 1 − cos 𝜃 − (𝑅2 − sen2 𝜃)2 ]} 𝑉𝑐
(36)
𝑅 =
𝑙
𝑎
(37)
em que 𝜃 é o ângulo do virabrequim; 𝑉(𝜃) é o volume do cilindro em função do ângulo
do virabrequim; 𝑅, apresentado na equação 37, é a razão entre o comprimento da biela e
o raio do virabrequim, como pode-se observar na figura 13; e 𝑉𝑐 é o volume morto no
interior do cilindro (volume quando o pistão se encontra no PMS), já comentado no item
3.1.
41
Figura 13: Simplificação gráfica de um cilindro
A figura 14 apresenta um gráfico do volume no interior do pistão em função do
ângulo do virabrequim calculado pela equação 36.
Figura 14: Volume em função do ângulo do virabrequim
PMS
PMI
𝑉 𝑑 𝑉
𝑙
𝑉 𝑐
𝑙 𝑠
𝑎
42
Utilizando os dados experimentais da pressão (P) e o valor calculado do volume
(V) no interior do cilindro, pode-se criar um gráfico de ln 𝑃 em função de ln 𝑉. A figura
15 apresenta um exemplo de uma curva de ln 𝑃 𝑥 ln 𝑉.
Figura 15: Gráfico de lnP x lnV.
O intervalo de compressão, antes do início da combustão, é uma função do
primeiro grau de inclinação 𝑘 = 𝑘𝑐𝑜𝑚𝑝, sendo que 𝑘 depende da composição dos gases
comprimidos. O mesmo ocorre para o intervalo de expansão, após o fim da combustão,
que também é uma função do primeiro grau de inclinação 𝑘 = 𝑘𝑒𝑥𝑝. 𝑘𝑒𝑥𝑝 é diferente de
𝑘𝑐𝑜𝑚𝑝 já que a composição dos gases antes e depois da combustão são diferentes.
Para calcular o início da combustão, o programa criado buscou identificar o
momento que a curva ln 𝑃 𝑥 ln 𝑉 deixa de ser uma equação do primeiro grau. O
procedimento utilizado foi o seguinte:
A- Utilizando um intervalo entre 10° após o fechamento da válvula de admissão
e 0° (Ponto Morto Superior) da curva ln 𝑃 𝑥 ln 𝑉, calcula-se o valor de R².
B- Se 𝑅2 ≥ 0,9999 , então o programa passa para a etapa seguinte. Caso
contrário, enquanto 𝑅2 < 0,9999, o valor de R² é recalculado, e a cada vez
que é recalculado o intervalo é diminuído. A diminuição do intervalo segue
43
o seguinte procedimento: enquanto mantêm-se a extremidade do intervalo
constante para 10° após o fechamento da válvula de admissão, a outra
extremidade diminui um ponto a cada interação (no sentido contrário à
compressão), levando à diminuição do intervalo analisado. O intervalo entre
dois pontos adjacentes é de 0,25°, já que este foi o intervalo de admissão dos
dados de pressão experimental.
C- Como o ajuste é feito utilizando um número muito grande de pontos, é
possível que um ou mais pontos após o início da combustão seja incluído no
intervalo encontrado no item B. Para garantir que isso não ocorra, após
encontrar o intervalo de pontos que permite um ajuste linear com 𝑅2 ≥
0,9999, este intervalo é diminuído de um valor arbitrário que foi escolhido
como sendo 10°.
D- Após a conclusão da etapa C, a etapa B é reaplicada mais uma vez e um novo
intervalo é encontrado.
E- Em seguida, o novo intervalo de pontos utilizado no item D para recalcular
o 𝑅2 passa por uma última análise. Para impedir que o ajuste linear seja feito
com poucos pontos, um critério é incluído. Caso o intervalo seja menor que
120 pontos, ou seja, 30°, o programa não salva os cálculos realizados até o
momento para o ciclo em questão e começa a calcular o próximo ciclo. Caso
contrário, os cálculos são salvos e o programa continua para as próximas
etapas.
F- É realizada uma regressão linear dos dados da pressão para o intervalo de
pontos calculado no item D. A figura 16 mostra um exemplo de um gráfico
da curva ln 𝑃 𝑥 ln 𝑉 sobreposta da curva gerada pela regressão linear.
44
Figura 16: Gráfico de lnP x lnV sobreposto da regressão linear na compressão
G- Em seguida é escolhido o intervalo de pontos entre 10° após o fechamento
da válvula de admissão e 40° após o PMS (ou seja, o intervalo começa na
região de apenas compressão e vai até um ponto muito depois do início da
combustão). Uma matriz é, então, criada para esse intervalo com a diferença
entre o valor de ln 𝑃 e a função do primeiro grau calculada em F, essa matriz
será chamada de 𝑣𝑎𝑟_𝑐𝑜𝑚𝑝.
H- Finalmente uma nova matriz será criada com a derivada da matriz
𝑣𝑎𝑟_𝑐𝑜𝑚𝑝, essa nova matriz será chamada de 𝑑𝑒𝑟_𝑣𝑎𝑟_𝑐𝑜𝑚𝑝.
I- Escolhendo apenas o intervalo definido em E, é encontrado o maior valor em
módulo de 𝑑𝑒𝑟_𝑣𝑎𝑟_𝑐𝑜𝑚𝑝. Esse ponto será chamado de
𝑚𝑎𝑥_𝑑𝑒𝑟_𝑣𝑎𝑟_𝑐𝑜𝑚𝑝. A figura 17 é um gráfico dos dados da matriz
𝑑𝑒𝑟_𝑣𝑎𝑟_𝑐𝑜𝑚𝑝 em função do ângulo do virabrequim. As duas linhas
pontilhadas representam o valor de ± 𝑚𝑎𝑥_𝑑𝑒𝑟_𝑣𝑎𝑟_𝑐𝑜𝑚𝑝. Como podemos
observar na figura, inicialmente há uma pequena variação dos dados em
torno do zero, ou seja, em torno da derivada zero. Em seguida há um
intervalo com pontos bastante distantes do zero, neste intervalo ocorre a
combustão.
45
Figura 17: Gráfico de 𝑑𝑒𝑟_𝑣𝑎𝑟_𝑐𝑜𝑚𝑝 em função do ângulo do virabrequim e das retas
de ± 𝑚𝑎𝑥_𝑑𝑒𝑟_𝑣𝑎𝑟_𝑐𝑜𝑚𝑝
J- O início da combustão foi definido como o primeiro ponto após o intervalo
calculado em E (no sentido crescente do valor do ângulo do virabrequim)
com valor da derivada de ln 𝑃 (calculada em H) 1,5 vezes maior que
𝑚𝑎𝑥_𝑑𝑒𝑟_𝑣𝑎𝑟_𝑐𝑜𝑚𝑝 . Caso ln 𝑃 > 1,5 ∙ 𝑚𝑎𝑥_𝑑𝑒𝑟_𝑣𝑎𝑟_𝑐𝑜𝑚𝑝 , mais um
teste é realizado: 40 pontos após o ponto em questão (no sentido da
compressão) também deverão ser 1,5 vezes maiores que
𝑚𝑎𝑥_𝑑𝑒𝑟_𝑣𝑎𝑟_𝑐𝑜𝑚𝑝.
As restrições descritas em J foram acrescentadas por dois motivos. Primeiro pois
a derivada de ln 𝑃 é muito elevada no início da combustão, ou seja, é bem maior que 1,5
vezes o erro dos pontos teóricos em relação à regressão linear de ln 𝑃 𝑥 ln 𝑉 no intervalo
de apenas compressão. O segundo motivo se dá devido ao fenômeno da injeção de
combustível. Quando tal fenômeno ocorre temos duas consequências: Primeiramente, é
criada uma turbulência maior no escoamento no interior do cilindro, e a segunda
consequência é que a composição do gás é alterada levemente (porém suficiente para que
haja um pequeno desvio da inclinação da regressão linear). Esse desvio ocorre pela
mudança do valor de 𝑘.
46
Foi observado empiricamente em diversas das curvas utilizadas que 1,5 vezes o
𝑚𝑎𝑥_𝑑𝑒𝑟_𝑣𝑎𝑟_𝑐𝑜𝑚𝑝 é suficiente para compensar as duas consequências geradas pela
injeção de combustível e garantir que o ponto do início da combustão seja corretamente
identificado.
4.5.2. Final da combustão
Como já mencionado, o fim da combustão é mais difícil de se identificar que o
seu início. Diferente do início, em que há uma liberação brusca de calor, o fim da
combustão é definido por uma diminuição gradual até o seu término. Por esse motivo dois
programas distintos foram criados e aplicados separadamente a fim de calcular o fim da
combustão. Em seguida os resultados foram comparados. Os dois métodos serão
comentados a seguir.
4.5.2.1. Método do erro
Os primeiros passos para calcular o fim da combustão são bastante similares ao
utilizado para calcular o seu início. Porém, alguns ajustes foram necessários. A seguir, o
procedimento de cálculo é detalhado. O programa que utiliza este método se encontra no
anexo I.
A- Utilizando o intervalo entre 10º antes da abertura da válvula de exaustão e 0º
(ponto morto superior) da curva ln 𝑃 𝑥 ln 𝑉, calcula-se o valor de R².
B- Analogamente ao procedimento do item 4.5.1 subitem B, se 𝑅2 ≥ 0,9998 (a
turbulência no processo de expansão é maior que na compressão, sendo
assim, foi necessário diminuir o critério do 𝑅2), então o programa passa para
a etapa seguinte. Caso contrário, enquanto 𝑅2 < 0,9998, o valor de R² é
recalculado, e a cada vez que é recalculado o intervalo é diminuído. A
diminuição do intervalo segue o seguinte procedimento: enquanto a
extremidade definida pelo ponto 10º antes da abertura da válvula de exaustão
se mantem constante, a extremidade anterior (sentido oposto à expansão) é
diminuída de um ponto a cada iteração, ou seja, a cada vez que 𝑅2 é
recalculado, o cálculo se inicia no ponto seguinte.
C- Análogo ao item C do item 4.5.1.
D- Análogo ao item D do item 4.5.1.
E- Análogo ao item E do item 4.5.1.
47
F- Análogo ao item F do item 4.5.1. Sendo que a figura 18 apresenta o gráfico
da curva ln 𝑃 𝑥 ln 𝑉 sobreposta da curva gerada pela regressão linear para a
exaustão.
Figura 18: Gráfico de lnP x lnV sobreposto da regressão linear na expansão
G- Em seguida, é escolhido o intervalo de pontos entre 20° antes do PMS e 10°
antes da abertura da válvula de exaustão (ou seja, o intervalo começa na
região anterior a combustão e se estende até um ponto que acredita-se ser
muito além dela). Uma matriz é então criada para esse intervalo com a
diferença entre o valor de ln 𝑃 e a função do primeiro grau calculada em F,
essa matriz será chamada de 𝑣𝑎𝑟_𝑒𝑥𝑝𝑎𝑛.
Diferente do cálculo do início da combustão, o final da combustão usará os dados
de 𝑣𝑎𝑟_𝑒𝑥𝑝𝑎𝑛 e não sua derivada.
H- Para o intervalo definido em D, é encontrado o valor mais distante da
regressão linear, esse valor será chamado de 𝑚𝑎𝑥_𝑣𝑎𝑟_𝑒𝑥𝑝𝑎𝑛. A figura 19
é um gráfico da matriz 𝑣𝑎𝑟_𝑒𝑥𝑝𝑎𝑛 e as linhas pontilhadas representam os
valores de ± 𝑚𝑎𝑥_𝑣𝑎𝑟_𝑒𝑥𝑝𝑎𝑛.
48
Figura 19: Gráfico de 𝑣𝑎𝑟_𝑒𝑥𝑝𝑎𝑛 em função do ângulo do virabrequim e das retas de
± 𝑚𝑎𝑥_𝑣𝑎𝑟_𝑒𝑥𝑝𝑎𝑛
I- O final da combustão é definido quando o ponto analisado tem seu respectivo
valor de 𝑣𝑎𝑟_𝑒𝑥𝑝𝑎𝑛 1,02 vezes maior que o módulo de 𝑚𝑎𝑥_𝑣𝑎𝑟_𝑒𝑥𝑝𝑎𝑛 e
quando 80 pontos (equivalente a um intervalo de 20º) antes do ponto em
questão (no sentido oposto à expansão) também apresentam seus respectivos
valores de 𝑣𝑎𝑟_𝑒𝑥𝑝𝑎𝑛 1,02 vezes maior que o módulo de 𝑚𝑎𝑥_𝑣𝑎𝑟_𝑒𝑥𝑝𝑎𝑛.
4.5.2.2. Método da menor derivada de ln 𝑃
O segundo método utilizado é proposto por VARAPRASADA et al. [48]. Este
método, determina que o fim da combustão pode ser encontrado pela menor derivada de
ln 𝑃 na etapa de Expansão. O programa baseado nesta teoria, que por sua vez se encontra
no anexo II, segue os passos abaixo.
A- Foi criada uma matriz da derivada de ln 𝑃 para o intervalo de 0° até 10° antes
do início da abertura da válvula de exaustão. Dentro desse intervalo foi
identificado o ponto com a menor derivada de ln 𝑃. Na figura 20 temos um
exemplo de um gráfico dessa matriz em função do ângulo do virabrequim.
49
Figura 20: Gráfico da derivada de lnP em função do ângulo do virabrequim
B- Como os dados experimentais são somados ao ruído proveniente da
admissão de dados, não é possível identificar claramente a derivada. Sendo
assim, um intervalo começando 7,5° antes da menor derivada e 7,5° depois
foi selecionado.
C- Foi realizada uma regressão polinomial do terceiro grau para o intervalo
selecionado em B. Na função gerada, foi identificado o seu mínimo e este foi
então considerado como o fim da combustão. A figura 21 mostra este
intervalo em uma curva para 100% de carga e o polinômio de grau três
gerado desses dados.
50
Figura 21: Gráfico da derivada de lnP em função do ângulo do virabrequim e de sua
regressão polinomial do terceiro grau.
4.5.3. Cálculo de IMEP
O IMEP (Pressão Média Efetiva Indicada), que é o trabalho líquido dividida pelo
volume deslocado (𝑉𝑑 ), considera o trabalho que é cedido pelo gás ao pistão, não
incluindo no cálculo as perdas mecânicas do motor. De acordo com HEYWOOD [17],
para calcular a potência indicada, primeiro deve-se calcular a área de B e de A, e em
seguida subtrair A de B. Essas áreas são representadas na figura 22. Vale apenas ressaltar
que a área de C é a intercessão de A e B.
A área A representa o trabalho realizado pelo gás sobre o pistão. Já a área B,
quando a pressão no interior do cilindro na etapa da admissão é menor que a da etapa de
exaustão, representa o trabalho do pistão sobre o gás, caso contrário representa o trabalho
realizado pelo gás sobre o pistão.
51
Figura 22: Parte de uma curva de pressão em função do ângulo do virabrequim
O cálculo de cada área é realizado pela equação 38 [17].
𝑊𝑖 = ∮ 𝑃 𝑑𝑉
(38)
O IMEP é então calculado pela equação 39:
𝐼𝑀𝐸𝑃 =
𝑊𝐴 − 𝑊𝐵
𝑉𝑑
(39)
Para calcular as áreas A e B foi utilizado a Regra do Trapézio. Essa regra calcula
a integral de uma função ligando cada dois pontos por uma função do primeiro grau, e
em seguida calculando a área abaixo desta reta para o intervalo entre os dois pontos. A
figura 23 apresenta um exemplo da integral de uma função sendo calculada pela Regra
do Trapézio [49].
Figura 23: Exemplo gráfico do cálculo da integral pela Regra do Trapézio
52
Uma área qualquer, 𝐴𝑞, calculada por este método utiliza a equação 40.
𝐴𝑞 =
1
2∙ (𝑓(𝑥𝑖−1) + 𝑓(𝑥𝑖)) ∙ (𝑥𝑖 − 𝑥𝑖−1)
(40)
onde o índice 𝑖 indica o i-ésimo ponto e 𝑓(𝑥𝑖) seu valor correspondente.
4.5.4. Cálculo de BMEP e FMEP
BMEP (Pressão Média Efetiva de Freio) é uma função da potência final
fornecida pelo motor (𝑃𝑏), após todas as perdas. Seu cálculo é realizado na equação 41
[17].
𝐵𝑀𝐸𝑃 =
𝑃𝑏 ∙ 𝑛𝑅
𝑉𝑑 ∙ 𝑁
(41)
Em que:
𝑛𝑅 = {1, 𝑝𝑎𝑟𝑎 𝑚𝑜𝑡𝑜𝑟𝑒𝑠 𝑑𝑒 2 𝑡𝑒𝑚𝑝𝑜;2, 𝑝𝑎𝑟𝑎 𝑚𝑜𝑡𝑜𝑟𝑒𝑠 𝑑𝑒 4 𝑡𝑒𝑚𝑝𝑜;
𝑁 é a rotação do motor; e
𝑉𝑑 é o volume deslocado.
Os valores de 𝑃𝑏 utilizados foram os valores da potência nominal do motor. As
equações de 42 a 45 apresentam as potências utilizadas para cada carga do motor.
𝑃𝑏 = 125 𝑘𝑊, 𝑝𝑎𝑟𝑎 25% 𝑑𝑎 𝑐𝑎𝑟𝑔𝑎 (42)
𝑃𝑏 = 250 𝑘𝑊, 𝑝𝑎𝑟𝑎 50% 𝑑𝑒 𝑐𝑎𝑟𝑔𝑎 (43)
𝑃𝑏 = 375 𝑘𝑊, 𝑝𝑎𝑟𝑎 75% 𝑑𝑒 𝑐𝑎𝑟𝑔𝑎 (44)
𝑃𝑏 = 500 𝑘𝑊, 𝑝𝑎𝑟𝑎 100% 𝑑𝑒 𝑐𝑎𝑟𝑔𝑎 (45)
FMEP (Pressão Média Efetiva de Atrito) é uma função da potência dissipada
pelo atrito dos componentes mecânicos do motor. Uma forma de calcular FMEP é
apresentada pela equação 46.
𝐹𝑀𝐸𝑃 = 𝐼𝑀𝐸𝑃 − 𝐵𝑀𝐸𝑃 (46)
53
4.5.5. 𝑃𝑚𝑎𝑥 e 𝜃𝑃𝑚𝑎𝑥
Os últimos cálculos realizados sobre as curvas experimentais foram os cálculos
da pressão máxima do ciclo (𝑃𝑚𝑎𝑥) e do ângulo desta pressão máxima (𝜃𝑃𝑚𝑎𝑥).
4.5.6. Tratamento dos dados
Os dados analisados neste trabalho foram submetidos ao Critério de Chauvenet
[50].
Supondo que diversas repetições de um teste realizados sobre as mesmas
condições são feitas, espera-se que os dados obedeçam uma distribuição Gaussiana.
Porém, algum fator externo pode influenciar em um ou mais resultados gerando dados
incompatíveis com o teste realizado. O Critério de Chauvenet elimina dados identificados
como suspeitos de terem sido influenciados.
De acordo com Chauvenet, espera-se que a probabilidade de um dado se
distanciar da média seja menor que 1 2𝑛⁄ , em que 𝑛 é o número total de dados [50]. Sendo
assim, os dados que estiverem fora deste intervalo serão considerados inválidos. A figura
24 exemplifica o que foi comentado até o momento, em que a área hachurada é a área
com os dados considerados validos por Chauvenet. Vale lembrar que os testes devem ser
realizados para as mesmas condições. Sendo assim, o Critério será aplicado
independentemente para as curvas de uma mesma carga.
Figura 24: Ilustração do critério de Chauvenet
54
A abscissa da figura acima é o desvio padrão em relação à média. Sendo assim,
para aplicar o Critério de Chauvenet, basta identificar o desvio padrão dos limites da área
hachurada, que por sua vez são iguais em módulo. O desvio padrão será chamado de
𝜎𝑚𝑎𝑥. Basta identificar os dados que estejam afastados da média por um desvio padrão
superior a 𝜎𝑚𝑎𝑥 e elimina-los. Para encontrar 𝜎𝑚𝑎𝑥 deve-se resolver a equação 47, sendo
que a equação 48 representa a Função Erro de Gauss (𝑒𝑟𝑓) [50].
1
2∗ erf (
𝜎𝑚𝑎𝑥
√2) =
1 − (1
2 ∗ 𝑛)
2
(47)
erf(𝜂1) =
1
√𝜋∫ 𝑒−𝜂2
+𝜂1
−𝜂1
𝑑𝜂 (48)
Os resultados calculados nos itens 4.5.1, 4.5.2.1 ou 4.5.2.2, 4.5.3 e 4.5.5 foram
submetidos independentemente ao Critério de Chauvenet, sendo um conjunto de dados
de uma mesma potência do motor por vez. Apenas os ciclos que tiveram os todos
resultados simultaneamente aprovados pelo Critério foram considerados validos para esse
trabalho.
Como é sugerido por Chauvenet, este critério foi aplicado apenas uma vez [50].
4.5.7. Curva de pressão representativa do ciclo
A fim de escolher uma curva de pressão experimental representativa do ciclo,
foram utilizados dois parâmetros da curva aplicando pesos diferenciados. Os parâmetros
escolhidos foram: diferença percentual da pressão máxima ( Δ%𝑃𝑚𝑎𝑥 ) e diferença
percentual do ângulo da pressão máxima ( Δ%𝜃𝑃𝑚𝑎𝑥) de cada ciclo. As diferenças
percentuais foram calculadas em relação às médias dos dados validados no item 4.5.6.
Foi então dado um peso duas vezes maior para Δ%𝑃𝑚𝑎𝑥 que para Δ%𝜃𝑃𝑚𝑎𝑥.
A pressão máxima recebeu um peso maior pois se 𝜃𝑃𝑚𝑎𝑥 for diferente do
calculado no item 4.5.6, este irá variar pelo menos 0,25º, que é o intervalo de admissão
dos dados experimentais, enquanto que a pressão pode apresentar qualquer valor. Sendo
assim, para que esse intervalo de 0,25º não apresente grande influência e restrinja demais
55
o número de curvas na escolha da curva representativa, foi dado essa diferença de peso.
A equação 49 apresenta o cálculo de 𝐴𝑈𝑋.
𝐴𝑈𝑋 =(2 ∗ |Δ%𝑃𝑚𝑎𝑥| + |Δ%𝜃𝑃𝑚𝑎𝑥
|)
3 (49)
Esta equação foi utilizada separadamente para cada curva experimental, a fim de
identificar uma curva de pressão representativa para cada carga no motor. A curva
escolhida foi a que apresentou menor valor de 𝐴𝑈𝑋 . A tabela 21 mostra o valor da
variável 𝐴𝑈𝑋 das curvas escolhidas para cada carga do motor.
Tabela 21: Valor da variável AUX para as curvas de pressão escolhidas de cada carga
Carga 𝑨𝑼𝑿 [%]
25% 0,09
50% 0,22
75% 0,11
100% 0,01
4.5.8. Resultados do programa das curvas de pressão
Na tabela 22 e 23 seguem os resultados das análises dos dois programas
descritos acima. São apresentados o início da combustão, o final da combustão, o IMEP,
o FMEP, a pressão máxima no cilindro (𝑃𝑚𝑎𝑥), número de curvas calculadas e número
de curvas validas pelo método de Chauvenet. Como já comentado, alguns valores citados
abaixo são uma média das curvas consideradas válidas pelo critério de Chauvenet, esses
valores estão acompanhados de seus respectivos desvios padrões (Desv. Pad.).
56
Tabela 22: Resultados do programa que inclui cálculo do fim da combustão pelo
Método do Erro
25% 50% 75% 100%
Valor
ou
Média
Desv.
Pad.
Valor
ou
Média
Desv.
Pad.
Valor
ou
Média
Desv.
Pad.
Valor
ou
Média
Desv.
Pad.
Curvas Calculadas 476 541 385 588
Curvas validas 469 536 383 573
Início da combustão [º] 0,02 0,25 -0,82 0,22 -1,05 0,20 -1,39 0,31
Final da combustão [º] 33,22 3,78 37,17 3,91 36,18 6,50 42,76 3,95
IMEP [bar] 7,19 0,12 12,04 0,12 17,07 0,17 21,81 0,10
BMEP [bar] 5,18 10,36 15,54 20,72
FMEP [bar] 2,01 1,68 1,52 1,08
𝑷𝒎𝒂𝒙 [bar] 72,9 0,8 101,2 1,3 135,0 1,8 160,6 1,6
Tabela 23: Resultados do programa que inclui cálculo do fim da combustão pelo
Método da menor derivada de lnP
25% 50% 75% 100%
Valor
ou
Média
Desv.
Pad.
Valor
ou
Média
Desv.
Pad.
Valor
ou
Média
Desv.
Pad.
Valor
ou
Média
Desv.
Pad.
Curvas Calculadas 425 0 380 0 582 0 586 0
Curvas validas 424 0 379 0 577 0 578 0
Início da combustão [º] -0,02 0,25 -0,82 0,20 -1,07 0,21 -1,40 0,31
Final da combustão [º] 30,77 0,98 32,41 1,05 35,67 0,62 38,44 0,52
IMEP [bar] 7,20 0,12 12,08 0,11 17,11 0,19 21,81 0,10
BMEP [bar] 5,18 10,36 0,00 15,54 20,72
FMEP [bar] 2,02 1,72 0,11 1,57 1,09
𝑷𝒎𝒂𝒙 [bar] 72,72 0,8 101,7 1,3 134,7 2,1 160,6 1,6
57
Como pode ser observado nas tabelas 22 e 23, há certa diferença entre o final da
combustão calculado com cada método. O motor foi simulado no BOOST alimentando a
equação de Wiebe com os dois dados do fim da combustão. Separadamente ambos os
casos foram ajustados assim como será descrito no item 5.1, ou seja, buscando igualar o
IMEP simulado ao experimental. Em seguida, ao comparar os valores das pressões
máximas simuladas e experimentais, pode-se observar que os casos que utilizaram os
dados calculados no item 4.5.2.1 apresentou valores mais próximos aos experimentais
que os calculados no item 4.5.2.2. Sendo assim, os dados da tabela 22 foram utilizados.
4.6. Cálculo da posição das válvulas de admissão e exaustão
Não foi encontrado no manual do motor informações sobre a abertura das
válvulas de admissão e exaustão. Sendo assim, foi necessário calcular esse valor.
Foi observado em um panfleto informativo do motor [51] de autoria do
fabricante, a empresa MAN, um desenho de corte lateral do pistão e mecanismos
auxiliares. A figura 25 apresenta esse desenho.
Figura 25: Figura de um corte do motor apresentando o mecanismo de abertura das
válvulas [51]
I
II
III IV
58
Em outros materiais informativos pode-se encontrar o valor nominal do diâmetro
da cabeça da válvula de admissão e de exaustão. O diâmetro da haste do corpo das
válvulas foi medido em duas peças sobressalentes disponíveis no LMT, sendo uma
válvula de admissão e uma de exaustão. Foi então feito a média de três medidas de cada
válvula. No manual pode ser encontrado o diâmetro do pistão do motor.
Como não foi especificando se os componentes do motor apresentados na figura
25 tinham suas dimensões proporcionais aos valores reais, nem qual seria essa proporção,
o procedimento a seguir foi realizado.
A figura foi importada para o programa AUTOCAD e as peças com as
geometrias encontradas acima foram medidas. Em seguida foram identificadas as
proporções entre cada geometria calculada no programa e as medidas em escala real. A
tabela 24 apresenta os dados geométricos, os dados calculados pelo programa
AUTOCAD e a razão entre esses dois valores. Será considerado que o deslocamento do
ponto D representa o deslocamento da válvula.
Tabela 24: Tabela com os valores reais de algumas peças, seu valor na figura
importada no AUTOCAD e a razão entre esses dois valores
Geometria
Valor em
escala real
[mm]
AUTOCAD 𝑨𝑼𝑻𝑶𝑪𝑨𝑫
𝑽𝒂𝒍𝒐𝒓 [𝒎𝒎]
Diâmetro da cabeça da
válvula de admissão 57 27,8279 0,488
Diâmetro da cabeça da
válvula de exaustão 54 26,3218 0,488
Diâmetro da haste das
duas válvulas 11,5 5,6004 0,487
Diâmetro do pistão 160 78,1367 0,487
Como pode-se observar pela quarta coluna da tabela a proporção entre os valores
apresentados pelo programa AUTOCAD e os valores de projeto/medidos, são
equivalentes, sendo assim, concluiu-se que as dimensões de todos os componentes
representados na figura seguirão a mesma proporção indicada pela quarta coluna. Será
utilizado a razão 0,487 calculada para o diâmetro do pistão como a razão entre os valores
medidos no AUTOCAD e os valores reais. Esta razão foi escolhida pois a dimensão do
59
diâmetro é bem maior que dos outros componentes, o que minimiza o erro da medida
realizada manualmente no AUTOCAD.
Como pode-se observar na figura 25, a abertura das válvulas é comandada pelo
seguintes mecanismos: um eixo com um excêntrico (I) que realiza uma volta a cada duas
do virabrequim; uma haste (II), que é apoiada sobre este eixo e realiza um movimento
vertical em função do contato do excêntrico; e um braço oscilante (balancim) (III), que
quando é deslocado em um sentido pela haste aciona a válvula (IV) impondo um
movimento no sentido contrário.
Uma simplificação do mecanismo de abertura da válvula é apresentada na figura
26, incluindo as dimensões cotadas no AUTOCAD para o momento da válvula fechada.
O objetivo é identificar quanto o ponto D se movimenta na vertical quando o ponto A é
deslocado devido à ação do excêntrico. O deslocamento vertical do ponto D será
considerado como o deslocamento das válvulas na mesma direção.
Figura 26: simplificação do mecanismo de abertura das válvulas e as medidas no
AUTOCAD para quando as válvulas se encontram fechadas
60
Para identificar o quanto o ponto A se movimenta, foi utilizada a figura 27. Esta
figura apresenta várias posições do ponto A equidistantes de um ângulo de 5°. Como o
mecanismo representado pelos círculos em vermelho não está necessariamente em
contato com o excêntrico na linha que conecta os dois centros, então foi realizado
manualmente a posição dos círculos em vermelho mantendo um intervalo de 5°, mas se
distanciando do excêntrico até que apenas um ponto do círculo vermelho estivesse em
contato com o excêntrico. Os valores de 𝑥 estão no anexo III, assim como o valor de 𝐷𝑥,
apresentado na figura 28, que é a diferença entre o valor de 𝑥 que está sendo analisado e
o raio do círculo primitivo do excêntrico (𝑅).
Figura 27: Diversas posições do excêntrico com intervalos de 5º
A figura 28 apresenta uma posição qualquer do mecanismo em um momento que
as válvulas se encontram abertas.
61
Figura 28: Posição do mecanismo de abertura das válvulas em um instante qualquer da
válvula aberta
Baseando-se nas figuras 26, 27 e 28 as equações a seguir foram criadas. A
equação 54 calcula o deslocamento do ponto D na vertical – 𝐷𝑒𝑠𝑙(𝐷). Para chegar na
equação 54 é necessário fazer algumas correlações geométricas do mecanismo de abertura
das válvulas. Essas correlações são apresentadas nas equações 50 a 53.
𝑅 = 162,87 ∗ sen 83 − 27,27 ∗ sen 2 + 27,27 ∗ sen 𝑇 − 𝐷𝑥 (50)
𝑆 = 19,14 − 27,27 ∗ (1 − cos 2) + 27,27 ∗ (1 − cos 𝑇) (51)
𝑄2 = 𝑅2 + 𝑆2 (52)
𝑃 = 30° − (𝑇 − 2°) (53)
𝐷𝑒𝑠𝑙(𝐷) = 43,77 ∗ (sen 30 − sen 𝑃) (54)
Para calcular o valor de 𝑇 foi utilizado o método de Newton. O método de
Newton é regido pela equação 55. Ela busca encontrar o valor de 𝑇 que satisfaça a
equação 𝑓(𝑇𝑛)=0 a partir de algumas iterações.
62
𝑇𝑛+1 = 𝑇𝑛 −𝑓(𝑇𝑛)
𝑓′(𝑇𝑛) (55)
onde 𝑛 indica a n-ésima iteração, 𝑓(𝑇𝑛) é uma função de 𝑇 de resultado zero e 𝑓′(𝑇𝑛) é
a derivada primeira de 𝑓(𝑇𝑛). No caso deste problema temos 𝑓(𝑇𝑛) regido pela equação
56.
𝑓(𝑇𝑛) = −𝑄2 + 𝑅2 + 𝑆2 = 0 (56)
Com auxílio de um programa criado no MATLab, que se encontra no anexo IV,
foi possível calcular o valor de T para diferentes valores de 𝐷𝑥 e, consequentemente,
calcular o 𝐷𝑒𝑠𝑙(𝐷), que é o quanto a válvula se desloca (abre). O critério de parada das
iterações foi: 𝑇𝑛+1 − 𝑇𝑛 ≤ 0,001°. Os dados de 𝐷𝑒𝑠𝑙(𝐷) estão no anexo III.
Há dois excêntricos na figura 29, sendo um responsável responsáveis pela
abertura das válvulas de admissão e o outro das de escape. O desenho dos dois excêntricos
em vermelho são iguais, apenas estão defasados de 92º. Pode-se observar que os dois
excêntricos desenhados se encaixam muito bem na figura base. Sendo assim, as aberturas
das válvulas em função do ângulo do virabrequim serão consideradas iguais, com uma
defasagem de 184º entre a válvula de exaustão e de admissão (já que cada volta do
excêntrico corresponde a duas voltas do virabrequim).
Figura 29: Imagem ampliada dos excêntricos de aberturas das válvulas de admissão e
exaustão, e uma geometria sobreposta criada no AUTOCAD.
63
5. Construção do modelo no BOOST
Para construção do modelo do motor no software BOOST são necessários
inicialmente os dados geométricos do motor. Muitos dados foram encontrados no manual
do motor [31], outros foram encontrados em outros materiais de autoria do fabricante [51,
52, 53, 54, 55] e os restantes foram aproximados a partir de medições realizadas no motor
instalado no LMT.
O manual apresenta informações de funcionamento do motor para quatro
condições de carga aplicada. Esses dados são apresentados na tabela 25.
Tabela 25: Dados de funcionamento do motor em função da carga
Parâmetros Carga
25% 50% 75% 100%
Temperatura na entrada da turbina [ºC] 320 410 450 500
Temperatura na saída da turbina [ºC] 275 325 320 320
Pressão na saída da turbina [bar] 1,002 1,002 1,002 1,002
Temperatura na entrada do compressor [ºC] 27 27 27 27
Pressão na entrada do compressor [bar] 1,0 1,0 1,0 1,0
Temperatura na saída do compressor [ºC] 70 120 152 184
Pressão na saída do compressor [bar] 1,4 2,1 2,7 3,5
Consumo específico de combustível [g/kW] 215 192 187 187
Consumo específico de ar [kg/kW] 11,2 8 7,5 6,8
A partir do esquema do motor descrito no manual e observado no LMT, o
modelo do motor foi construído no BOOST e é apresentado na figura 30.
64
Figura 30: Imagem da tela do programa BOOST
5.1. Ajustes
Alguns parâmetros do BOOST devem ser ajustados para adequar o modelo
simulado aos dados experimentais.
5.1.1. Turbo compressor
O turbo compressor apresenta alguns parâmetros que podem ser ajustados pelo
usuário. São eles: eficiência do compressor, eficiência da turbina, eficiência mecânica do
turbo compressor e correção do fluxo de massa. Os parâmetros foram ajustados
iterativamente buscando equiparar os dados experimentais da vazão mássica de ar,
temperatura e pressão do ar comprimido. Como não foi possível encontrar um ajuste que
equiparasse os dados simulados aos experimentais foi dada uma importância maior à
vazão mássica de ar comprimido. A tabela 26 apresenta uma comparação entre os dados
simulados e experimentais do gás na saída do compressor após o ajuste dos parâmetros
do turbo compressor.
65
Tabela 26: Informações experimentais e simuladas do motor
Carga Pressão [bar] Temperatura [K] Vazão mássica [kg/s]
SIM EXP DIF% SIM EXP DIF% SIM EXP DIF%
25% 1,51 1,40 7,86 340 343 -0,87 0,3866 0,3889 -0,6
50% 2,18 2,10 3,81 393 393 0,00 0,5617 0,5556 1,1
75% 2,97 2,70 10,00 414 393 5,34 0,7495 0,7813 -4,1
100% 3,82 3,50 9,14 495 458 8,08 0,9456 0,9444 0,13
As siglas SIM, EXP e DIF% se referem respectivamente à dados simulados,
dados experimentais e diferença percentual entre SIM e EXP.
5.1.2. IMEP
Com o objetivo de equiparar o IMEP experimental ao IMEP simulado, a
transferência de calor do gás no interior do cilindro com a parede foi ajustada. Dois
parâmetros que podem ser ajustados afetam a taxa de transferência de calor do gás com a
parede do cilindro: as temperaturas da parede e o fator multiplicado do coeficiente de
transferência de calor ℎ𝐴𝑉𝐿 (FMCTC). Foram determinadas temperaturas fixa para a parede
do cilindro, sendo que o valor dessas temperaturas foram baseados no livro de
HEYWOOD [17]. FMCTC foi ajustado de forma iterativa para cada carga do motor. Tanto
o valor das temperaturas fixadas quanto de FMCTC são apresentadas na tabela 27.
Tabela 27: Dados usados para ajustar o modelo de transferência de calor
Carga Temperatura da parede
do cilindro no PMI
Temperatura da parede das
outras partes do cilindro FMCTC
25% 90 130 0,58
50% 90 130 0,81
75% 90 130 1,50
100% 90 130 1,30
66
Como não foi alterada a temperatura da parede, que na realidade deve mudar
com a carga do motor, o fator FMCTC também terá o papel de compensar essa diferença.
A tabela 28 apresenta a comparação entre os dados simulados e os experimentais
do IMEP, quando utilizados os valores dos parâmetros apresentados na tabela 27 e diesel
puro como combustível.
Tabela 28: comparação do IMEP simulado e experimental após ajuste do modelo do
coeficiente de calor
Carga IMEP [bar]
EXP SIM DIF%
25% 7,19 7,16 0,42
50% 12,04 12,07 -0,30
75% 17,07 16,71 2,11
100% 21,81 21,83 -0,09
5.1.3. FMEP
Por último, o FMEP foi ajustado. Para conciliar valor de FMEP simulado com o
experimental, há apenas um parâmetro a ser ajustado: o parâmetro multiplicativo do
FMEP (MFMEP). Este parâmetro foi ajustado de forma iterativa em função da carga
aplicada no motor. Os valores de FMEP simulado (SIM), experimental (EXP) e a
diferença percentual entre ambos (DIF%) são apresentados na tabela 29.
Tabela 29: Valores de FMEP experimentais e simulados após ajuste de MFMEP
Carga FMEP [bar]
MFMEP EXP SIM DIF%
25% 0,790 2,01 2,02 -0,5
50% 0,538 1,68 1,72 -2,4
75% 0,405 1,52 1,54 -1,3
100% 0,242 1,08 1,09 -0,9
67
Foi realizado um ajuste linear dos valores do parâmetro MFMEP em função da
carga aplicada no motor. Como pode-se observar na figura 31, há uma forte correlação
linear desse parâmetro em função da carga aplicada no motor
Figura 31: Dados de MFMEP e ajuste linear desses dados
O aumento da temperatura no interior do cilindro gera uma maior dilatação do
pistão, o que causa uma maior pressão entre os anéis de segmento e a parede do pistão.
Este fenômeno, individualmente, seria responsável por aumentar o atrito entre os anéis
de segmento e o pistão, consequentemente aumentando o FMEP. Porém, MACEK et al.
[56] e TOMANIK et al. [57], mostraram que enquanto o óleo se encontra na região de
lubrificação hidrodinâmica [58], o aumento da temperatura causa queda do coeficiente de
atrito entre o óleo e a as partes metálicas em contato com ele. Isto ocorre devido à
diminuição da viscosidade, e consequente redução da espessura de óleo. Sendo assim,
FMEP diminui.
R² = 0,9808
0
0,1
0,2
0,3
0,4
0,5
0,6
0,7
0,8
0,9
20 40 60 80 100
MFM
EP
% da carga do motor
M_FMEP
Ajuste Linear
68
6. Resultados
6.1. Validação da simulação
Os dados experimentais foram gerados para o motor trabalhando a 25%, 50%,
75% e 100% de carga nominal com o uso de diesel marítimo. Diferente do óleo diesel
automotivo vendido hoje no Brasil, o diesel marítimo não apresenta biodiesel, e por esse
motivo ele será representado por E0B0. Na tabela 30 são apresentados os valores
experimentais (EXP) simulados (SIM) e a suas diferenças percentuais (%Dif) para os
seguintes dados: IMEP, vazão mássica de ar que passa pelo compressor (𝑉𝑎𝑟) e a pressão
máxima atingida dentro da câmara de combustão (𝑃𝑚𝑎𝑥). Os valores experimentais são
referentes às curvas características de cada carga, com exceção de 𝑉𝑎𝑟 que foi calculado
a partir de dados contidos no manual.
Tabela 30: Valores simulados e experimentais do diesel puro (E0B0)
Carga IMEP (bar) 𝑽𝒂𝒓 (kg/s) 𝑷𝒎𝒂𝒙 (bar)
EXP SIM %Dif EXP SIM %Dif EXP SIM %Dif
25% 7,19 7,16 0,42 0,3889 0,3866 0,59 73,5 76,8 -4,49
50% 12,04 12,076 -0,30 0,5556 0,5617 -1,10 101,2 96,7 4,45
75% 17,07 16,71 2,11 0,7812 0,7496 4,05 136,4 143,4 -5,13
100% 21,81 21,83 -0,09 0,9444 0,9456 -0,13 160,2 163,3 -1,94
As figuras 32, 33, 34 e 35 apresentam gráficos para as quatro cargas aplicadas
no motor. Nos gráficos são comparadas as curvas de combustão simuladas (SIM) e
experimentais (EXP). Vale lembrar que as curva de pressões experimentais foram
encontradas no item 4.5.7.
69
Figura 32: Curvas de pressão simulada e experimental para 25% de carga do motor
Figura 33: Curvas de pressão simulada e experimental para 50% de carga do motor
0
10
20
30
40
50
60
70
80
90
-150,00 -100,00 -50,00 0,00 50,00 100,00
25% de carga
SIM EXP
0
20
40
60
80
100
120
-150,00 -100,00 -50,00 0,00 50,00 100,00
50% de carga
SIM EXP
70
Figura 34: Curvas de pressão simulada e experimental para 75% de carga do motor
Figura 35: Curvas de pressão simulada e experimental para 100% de carga do motor
Pode-se observar pelas figuras apresentadas que o modelo simulado apresenta
resultados muito próximos do experimental. Sendo assim, o modelo simulado é válido
pelo menos para as condições de operação descritas.
A diferença das curvas no estágio final da compressão e início da combustão,
em que a curva experimental apresenta valores da pressão maiores que a simulada,
0
20
40
60
80
100
120
140
160
-120,00 -70,00 -20,00 30,00 80,00
75% de carga
SIM EXP
0
20
40
60
80
100
120
140
160
180
-150,00 -100,00 -50,00 0,00 50,00 100,00
100% de carga
SIM EXP
71
também ocorreu em outros trabalhos que utilizaram o software BOOST. MELO [5]
encontrou esta característica da curva simulada ao utilizar o modelo de combustão Wiebe
duas zonas para um motor Flex-Fuel. MATIAS et al. [59] também encontraram esta
característica da curva ao simular o mesmo motor utilizando o modelo fractal de
combustão, porém ao comparar com o trabalho de MELO pode-se observar uma
diminuição deste fenômeno em relação às curvas experimentais.
6.2. Consumo de combustível para diferentes misturas
6.2.1. 100% de Carga nominal do motor
Os parâmetros 𝑎 e 𝑚 de Wiebe buscam computar diferentes condições de
operação do motor. Neste trabalho, 𝑎 foi mantido constante a fim de determinar o
percentual de combustível queimado pelo modelo de combustão de Wiebe. O parâmetro
𝑚 foi aproximado como sendo dependente apenas do combustível uma vez que foi
definido em 4.3.1. Pode-se observar neste item que a aproximação da curva da fração de
massa queimada para o modelo de Wiebe duas zonas não apresentou grandes mudanças
do parâmetro 𝑚, o que reforça a suposição que o modelo validado para o combustível
E0B0 deve ser válido para os outros combustíveis testados.
Mantendo todos os parâmetros ajustados para a simulação de E0B0 da vazão de
combustível injetada no cilindro, simulações para as misturas E1B7, E1B10, E1B20,
E2B7, E2B10, E2B20, E5B7, E5B10 e E5B20 foram realizadas. A vazão de combustível
foi alterada manualmente até que fosse encontrado o IMEP da curva simulada de 100%
de carga do motor. Na tabela 31 encontra-se aos resultados do IMEP, vazão mássica de
combustível (𝑉𝑐𝑜𝑚𝑏) e vazão mássica de ar no compressor (𝑉𝑎𝑟) para E0B0 e cada
mistura. %Dif indica a diferença percentual do valor de cada parâmetro em relação aos
valores encontrados pela simulação de E0B0.
72
Tabela 31: Valores simulados com todas as misturas de combustíveis e comparação
com os valores simulados para E0B0
IMEP
[bar] %Dif 𝑽𝒂𝒓 [
𝒌𝒈
𝒔] %Dif 𝑽𝒄𝒐𝒎𝒃 [
𝒌𝒈
𝒔] %Dif
E0B0 21,827 - 0,94561 - 0,02597 -
E1B7 21,826 0 0,9446 -0,11 0,02621 0,9
E1B10 21,829 0,01 0,94451 -0,12 0,02626 1,1
E1B20 21,824 -0,01 0,94452 -0,12 0,02642 1,7
E2B7 21,824 -0,01 0,94442 -0,13 0,02631 1,3
E2B10 21,8239 -0,01 0,9446 -0,11 0,02636 1,48
E2B20 21,825 -0,01 0,94438 -0,13 0,02653 2,11
E5B7 21,823 -0,02 0,9446 -0,11 0,02663 2,48
E5B10 21,8224 -0,02 0,94453 -0,11 0,02669 2,68
E5B20 21,826 0 0,94446 -0,12 0,02687 3,33
A figura 36 apresenta uma representação gráfica considerando os percentuais de
biodiesel, etanol e aumento do consumo de combustível em relação a E0B0.
Figura 36: Representação gráfica do aumento percentual de consumo de combustível
em função do percentual de etano e de biodiesel da mistura
7
10
20
0,00
0,50
1,00
1,50
2,00
2,50
3,00
3,50
1
2
5
Teor de biodiesel [%]
% d
e au
men
to d
e V
_co
mb
Teor de ETANOL [%]
0,00-0,50 0,50-1,00 1,00-1,50 1,50-2,00 2,00-2,50 2,50-3,00 3,00-3,50
73
Como pode ser observado na tabela 31 e na figura 36 que tanto o aumento de
etanol quanto o de biodiesel resultam em um maior consumo de combustível. Outro ponto
importante que deve ser destacado é que o aumento de etanol apresenta um peso muito
maior no aumento do consumo de combustível que o biodiesel. Este comportamento das
misturas está diretamente relacionado com o PCI (Poder Calorífico Inferior) de cada
substância pura. Como pode ser observado nos itens 4.2.1, 4.2.2 e 4.2.3 o Poder Calorífico
do biodiesel e do etanol são inferiores ao do diesel, porém o do biodiesel é muito mais
próximo do diesel que o do etanol. Portanto, será necessário injetar uma massa de
combustível maior para as misturas do que para o diesel puro (E0B0), para não haver
perda de potência no motor
6.2.2. 50% de carga
Assim como foi realizado para no item 6.2.1, foi simulado o motor para 50% de
carga, variando a vazão de combustível para as diferentes misturas, a fim de encontrar o
IMEP equivalente à curva simulada para E0B0. Na tabela 32 é apresentado os resultados
do IMEP, vazão mássica de combustível (𝑉𝑐𝑜𝑚𝑏) e vazão mássica de ar no compressor
(𝑉𝑎𝑟) para E0B0 simulada e para cada mistura. %Dif indica a diferença percentual de
cada parâmetro em relação aos dados simulados de E0B0. Os combustíveis simulados
foram: E2B7, E2B20, E5B7, e E5B20
Tabela 32 Valores simulados de algumas misturas de combustíveis e comparação com
os valores simulados para E0B0
IMEP
[bar] %Dif 𝑽𝒂𝒓 [
𝒌𝒈
𝒔] %Dif 𝑽𝒄𝒐𝒎𝒃 [
𝒌𝒈
𝒔] %Dif
E0B0 12,071 - 0,5617 - 0,01334 -
E2B7 12,07 0,01 0,5617 0 0,0135 1,25
E2B20 12,069 -0,02 0,5604 -0,23 0,01362 2,11
E5B7 12,063 -0,07 0,5615 -0,03 0,01366 2,42
E5B20 12,071 0 0,5687 1,25 0,01379 3,39
Analogamente ao que foi observado no item 6.2.1, a vazão de combustível
aumentou com o aumento no teor de etanol e de biodiesel. Sendo que a influência do PCI
74
apresentou correlação com o aumento do consumo de combustível similar ao que foi
discutida em 6.2.1.
Observando o aumento percentual de consumo de combustível nas tabelas 31
(100% de carga) e 32 (50% de carga) entre os mesmos combustíveis, encontra-se que a
maior diferença de %Dif foi de 0,06%. Pode-se concluir então que o consumo de
combustível não apresenta dependência relevante com relação a carga aplicada ao motor.
75
7. Conclusões
Este trabalho apresentou o desenvolvimento de um modelo computacional de
um motor marítimo do tipo diesel. Foi desenvolvido no software BOOST da empresa
AVL utilizando o modelo de combustão Wiebe duas zonas e o modelo de transferência
de calor desenvolvido pela empresa AVL. O modelo foi validado com dados
experimentais para 25%, 50%, 75% e 100% da carga nominal utilizando óleo diesel
marítimo.
Ajustando os parâmetros do software, a fim de encontrar o IMEP calculado dos
dados experimentais, foi verificada pequena variação quando comparada as curvas de
pressão simuladas e experimentais. Sendo assim, foi possível validar o modelo
computacional para as quatro condições de operação do motor citadas.
Foram simuladas diferentes misturas de diesel adicionadas de biodiesel e etanol.
O enfoque deste estudo foi identificar a diferença do consumo de combustível das
diferentes misturas compostas, procurando-se manter a potência do motor. Foi
identificado então aumento no consumo de todas as misturas em relação ao combustível
base, diesel marítimo (sem adição de etanol nem biodiesel). Observou-se que o etanol
tem uma influência muito maior que o biodiesel em aumentar o consumo, e esse
fenômeno foi atribuído ao valor do PCI do etanol que é consideravelmente inferior ao do
óleo diesel. Analogamente o biodiesel também apresenta um PCI menor que o do diesel,
porém bem mais elevado que o etanol, por isso seu impacto no aumento do consumo não
é tão grande quanto o do etanol.
O resultado do consumo de combustível de algumas misturas foi comparado para
100% e 50% de carga. Concluiu-se que a condição de operação do motor apresentou
mínima influência no consumo, levando a concluir que o principal fator responsável pelo
aumento do consumo é a composição das misturas.
Sugestões para projetos futuros:
Simulação da combustão considerando o atraso de ignição como função
dos combustíveis;
Simulação utilizando o modelo de combustão Double Wiebe;
Investigar emissões em função das misturas;
Validação das curvas de pressão simuladas para as diferentes misturas;
76
8. Referências Bibliográficas
[1] http://www.mwm.com.br/site.aspx/Newsletters-Detalhe?codNoticia=
TZxgBRp/Zkg= , acessado em 27/06/2014 às 19:15 horas;
[2] EUROPEAN RENEWABLE ENERGY COUNCIL, “Renewable Energy Technology
Roadmap 20% by 2020”, 2008.
[3] AMEZAGA, J.M., BOYES, S.L. AND HARRISON, J.A., “Biofuels Policy in the
European Union”, 7th International Biofuels Conference, New Delhi, 2010
[4] http://www.biofuelstp.eu/global_overview.html, acessado em 27/06/2014 às 19:37
horas;
[5] MELO, T. C. C., Análise Experimental e Simulação Computacional de um Motor
Flex Operando com Diferentes Misturas de Etanol Hidratado na Gasolina, Tese de D.Sc.,
COPPE/UFRJ, Rio de Janeiro, RJ, Brasil, 2012.
[6] ENVIRONMENTAL PROTECTION AGENCY (EPA), “A Comprehensive Analysis
of Biodiesel Impacts on Exhaust Emissions”, Draft Technical Report, United States,
EPA420-P-02-001, 2002.
[7] ROCHA, H. L., “Estudo Experimental do Uso de Misturas Combustíveis Ternárias
em um Motor CFR/ASTMCETANO, Projeto Final de Graduação, Universidade Federal
do Rio de Janeiro, Rio de Janeiro, RJ, Brasil, 2014.
[8] SHI, X., PANG, X., UM, Y., HE, H., SHUAI, S., WANG, J., CHEN, H., LI, R.,
“Emission reduction potential of using ethanol–biodiesel–diesel fuel blend on a heavy-
duty diesel engine”, Atmospheric Environment, 2005.
[9] AL-HASSAN, M., MUJAFET, H., AL-SHANNAG, M., “An Experimental Study on
the Solubility of a Diesel-Ethanol Blend and on the Performance of a Diesel Engine
Fueled with Diesel-Biodiesel - Ethanol Blends”, Jordan Journal of Mechanical and
Industrial Engineering, v. 6, n. 2, pp.147 - 153 , 2012.
77
[10] EUROPEAN PARLIAMENT, EUROPEAN COUNCIL, “Directive 2009/28/EC”,
Official Journal of the European Union L 140, pp. 16-47, 2009.
[11] KARAVALIS, G., FONTARAS, G., BAKEAS, E., STOURNAS, S., “Effect of
Biodiesel Origin on the Regulated and PAH Emissions from a Modern Passenger Car”,
SAE Technical Paper, 2011.
[12] ROBBINS, C., HOEKMAN, S. K., CENICEROS, E., NATARAJAN, M., “Effects
of Biodiesel Fuels Upon Criteria Emissions”, SAE Technical Paper 2011-01-1943, 2011
[13] ACYOLY, A. F., FERNANDES, M. R., NETO, A. A. D., et al., “Investigação de
Propriedades do Sistema Etanol-Diesel-Biodiesel Aditivado com Tensoativo para
Aplicação como Combustível Alternativo”, 5° Congresso Brasileiro de Pesquisa e
Desenvolvimento em Petróleo e Gás, Fortaleza, Ceará, 15-22 março 2009.
[14] AGÊNCIA NACIONAL DO PETRÓLEO, GÁS NATURAL E
BIOCOMBUSTÍVEIS, “Resolução ANP Nº 7, DE 19.3.2008 - DOU 20.3.2008”, 2008.
[15] AGÊNCIA NACIONAL DO PETRÓLEO, GÁS NATURAL E
BIOCOMBUSTIVEIS, Resolução ANP Nº 50, DE 23.12.2013 - DOU 24.12.2013, 2013.
[16] LOGAHATHAN, S., TAMILPORAI, P., VIJAYAN, K., et al, “Simulation and
Analysis of Effect of Oxygenate Blended Diesel on Combustion and Performance in
Turbocharged Diesel Engine”, SAE Technical Paper 2007-01-2019, 2007.
[17] HEYWOOD, J. B., Internal Combustion Engine Fundamentals. 1 ed. Nova Deli,
Tata McGraw-Hill, 1988.
[18] ANNAND W. J. D.,”Heat transfer in the cylinders of reciprocating IC engines”
Proc.I.Mech.E. London, Vol.177, No.36 (1963).
[19] WOSCHNI, G., “A universally applicable equation for the instantaneous heat
transfer coefficient in the internal combustion engine”, SAE Technical Papers, 670931,
1967.
78
[20] MACHADO, G. B., Metodologias para Desenvolvimento de Combustíveis e
Determinação da Velocidade de Propagação da Chama em Motores de Ignição por
Centelha, Tese de D.Sc., PUC-Rio, Rio de Janeiro, RJ, Brasil, 2012
[21] STONE, R. “Introduction to internal combustion engines”, Warrendale: SAE, 1995
[22] BARROS, J. E. M., “Estudo de motores de combustão interna aplicando linguagem
orientada a objetos” Doutorado (D.Sc.) – Departamento de Engenharia Mecãnica,
Universidade Federal de Minas Gerais, Belo Horizonte, 2003
[23] POULOS S.G., HEYWOOD G.B., “The Effect of Chamber Geometry on Spark-
Ignition Engine Combustion”, SAE Paper 830334, 1983
[24] MELO, T. C.C., Modelagem Termodinâmica de um Motor Otto tipo Flex-Fuel,
Funcionando com Gasolina, álcool e Gás Natural, Tese de M.Sc., COPPE/UFRJ, Rio de
Janeiro, RJ, Brasil, 2007
[25] LANZAFAME, R., MESSINA, M., 2003, “ICE gross heat release strongly
influenced by specific heat ratio values”, International Journal of Automotive
Technology, Vol. 4, N.3, pp. 125-133.
[26] http://www.britannica.com/EBchecked/topic/434878/Nikolaus-August-Otto,
acessado em 03/04/2014 às 22:07 horas
[27] http://www.deutsches-museum.de/en/collections/machines/power-
engines/combustion-engines/diesel-engines/the-first-diesel-engine-1897/, acessado em
03/04/2014 às 21:23 horas
[28] AMERICAN SOCIETY FOR TESTING AND MATERIALS. ASTM D2699:
Standard Test Method for Research Octane Number of Spark-Ignition Engine Fuel.
Estados Unidos da América, 2013.
79
[29] AMERICAN SOCIETY FOR TESTING AND MATERIALS. ASTM D2700:
Standard Test Method for Motor Octane Number of Spark-Ignition Engine Fuel. Estados
Unidos da América, 2013.
[30] AMERICAN SOCIETY FOR TESTING AND MATERIALS. ASTM D613:
Standard Test Method for Cetane Number of Diesel Fuel Oil. Estados Unidos da América,
2013.
[31] L16/24 Project Guide, MAN B&W Diesel.
[32] http://www.planeta.coppe.ufrj.br/artigo.php?artigo=1183, acessado em 27/06/2014
às 01:09 horas
[33] JANAF Thermodynamic tables, 2d es, NSRDS-NB537, U.S. National Bureau of
Standards, June 1917.
[34] http://webbook.nist.gov/cgi/cbook.cgi?ID=C7440440&Units=SI, acessado em
01/04/2014 às 18:01 horas
[35] http://webbook.nist.gov/cgi/cbook.cgi?ID=C12385136&Units=SI, acessado em
01/04/2014 às 18:02 horas
[36] http://webbook.nist.gov/cgi/cbook.cgi?Formula=O&NoIon=on&Units=SI,
acessado em 01/04/2014 às 18:02 horas
[37] FILHO, F. F., Desempenho de um Motor Diesel Operando no Modo Bi-Combustível
Diesel-Álcool, Departamento de Engenharia Mecânica, PUC-Rio, Rio de Janeiro, RJ,
Brasil, 2008.
[38] Resolução ANP Nº 7, DE 9.2.2011 - DOU 10.2.2011 – RETIFICADA DOU
14.4.2011, Agência Nacional do Petróleo, Gás Natural e Biocombustíveis, 2011.
[39] BORGNAKKE, C., SONNTAG, R .E., Fundamentos as Termodinâmica. 7 ed. São
Paulo, Blucher, 2009
80
[40] VLASSOV, D. VIII Curso de Especialização em Motores e Combustíveis, 2008.
[41] LIN, CY., LIN, HA., “Diesel engine performance and emission characteristics of
biodiesel produced by the peroxidation process”, Fuel, v.85, n.3, pp.298-305, 2006
[42] LIU, HP., STRANK, S., “Emissions and Testing oh Neat Biodiesel Fuels”, SAE
Technical Paper 2010-90038, 2010
[43] http://www2.camara.leg.br/camaranoticias/noticias/CONSUMIDOR/469676-MP-
647-AUMENTA-PARA-7-O-PERCENTUAL-DE-BIODIESEL-NO-OLEO-
DIESEL.html --- 16-06-14, acessado em 03/04/2014 às 21:07 horas.
[44] http://www.biodiesel.org/using-biodiesel/oem-information/oem-statement-
summary-chart, acessado em 03/04/2014 às 21:20 horas
[45] GOGOI, T. K., SARMA, A. K., MISRA, P. S., SYED, T. H., “Combustion Analysis
of Jatropha Methyl Ester and Its Ethanol and Acetone Blends In a Diesel Engine”,
International Journal of Emerging Technology and Advanced Engineering, Vol. 3,
Special Issue 3: ICERTSD 2013, pages 51-57, Feb 2013
[46] WOSCHNI, G., HUBER, K., "The Influence of Soot Deposits on Combustion
Chamber Walls on Heat Losses in Diesel Engines," SAE Technical Paper 910297, 1991
[47] AVL, “AVL BOOST software version 2013 – Theory”, Graz, Austria, 2013
[48] VARAPRASADA R., NEHRU K. V. K., GANESAN V., “Evaluation of S.I. engine
combustion parameters: a new approach”, Combustion Science and Technology, v.89,
Issue 1-4, pp. 47-55, 1993.
[49] RUGGIERO M., LOPES V., Cálculo Numérico, Aspectos Teóricos e
Computacionais, 2 Ed., São Paulo, Pearson, 1988.
81
[50] HOLMAN, J. P., Experimental Methods for Engines. 6 ed. United States, McGraw-
Hill, 1966.
[51] Service Experience, Small Bore Four-stroke Engines, 5510-0118-00ppr Aug 2012,
MAN B&W Diesel, Dinamarca, 2012.
[52] L16/24 Project Guide – Marine, Complete manual, MAN B&W Diesel, 2013.
[53] Sedef 146, Operating Results, MAN B&W Diesel, 2007.
[54] Service Letter SL09-510/SBJ, MAN B&W Diesel, Abril 2009.
[55] MAN L16/24, GenSet power for marine applications, 1510-0106-02ppr June 2010,
MAN B&W Diesel, 2010.
[56] MACEK, J., FUENTE, D., EMRICH, M., "A Simple Physical Model of IC E
Mechanical Losses," SAE Technical Paper 2011-01-0610, 2011.
[57] TOMANIK, E., FERRARESE, ANDRE., “Low Friction Ring Pack for Gasoline
Engines”, ASME technical 2006-1566, 2006.
[58] BUDYNAS, R., NISBETT, J., Elementos de Máquinas de Shigley, 8 ed. United
States, McGraw-Hill, 2011.
[59] MATIAS, F. A. S., MELO T. C. C., MACHADO, G. B., “Comparação entre os
Modelos de Combustão Wiebe e Fractal para Simulação Computacional de um Motor
Flex Operando com Diferentes Misturas de Gasolina-Etanol”, SIMEA, 2014.
82
Anexo I: Programa para calcular dados da curva de pressão e utilizando o método
descrito no item 4.5.2.1.
%% Limpando todas as variaveis e importando os dados clear;
% Importando as curvas de pressão de três tomadas diferentes e em
seguida juntando as curvas de mesma carga mas de tomadas diferentes em
uma única matriz
P25_1=importdata('Tomada 01\P25.txt'); P25_1=P25_1.data; P25_2=importdata('Tomada 02\P25.txt'); P25_2=P25_2.data; P25_3=importdata('Tomada 03\P25.txt'); P25_3=P25_3.data;
P50_1=importdata('Tomada 01\P50.txt'); P50_1=P50_1.data; P50_2=importdata('Tomada 02\P50.txt'); P50_2=P50_2.data; P50_3=importdata('Tomada 03\P50.txt'); P50_3=P50_3.data;
P75_1=importdata('Tomada 01\P75.txt'); P75_1=P75_1.data; P75_2=importdata('Tomada 02\P75.txt'); P75_2=P75_2.data; P75_3=importdata('Tomada 03\P75.txt'); P75_3=P75_3.data;
P100_1=importdata('Tomada 01\P100.txt'); P100_1=P100_1.data; P100_2=importdata('Tomada 02\P100.txt'); P100_2=P100_2.data; P100_3=importdata('Tomada 03\P100.txt'); P100_3=P100_3.data;
for y=1:3 if y==1; P25=P25_1; else if y==2 for x=2:size(P25_2,2) P25(:,size(P25,2)+1)=P25_2(:,x); end; else for x=2:size(P25_3,2) P25(:,size(P25,2)+1)=P25_3(:,x); end; end end
if y==1; P50=P50_1; else if y==2 for x=2:size(P50_2,2)
83
P50(:,size(P50,2)+1)=P50_2(:,x); end; else for x=2:size(P50_3,2) P50(:,size(P50,2)+1)=P50_3(:,x); end; end end
if y==1; P75=P75_1; else if y==2 for x=2:size(P75_2,2) P75(:,size(P75,2)+1)=P75_2(:,x); end; else for x=2:size(P75_3,2) P75(:,size(P75,2)+1)=P75_3(:,x); end; end end
if y==1; P100=P100_1; else if y==2 for x=2:size(P100_2,2) P100(:,size(P100,2)+1)=P100_2(:,x); end; else for x=2:size(P100_3,2) P100(:,size(P100,2)+1)=P100_3(:,x); end; end end
end;
%% Loop para calcular as curvas de pressão de diferentes cargas for carga=1:4 q=0;
if carga==1 size_xlsx=size(P25,2); P_IMEP=P25; BMEP=500*0.25; end; if carga==2 size_xlsx=size(P50,2); P_IMEP=P50; BMEP=500*0.50; end; if carga==3 size_xlsx=size(P75,2); P_IMEP=P75; BMEP=500*0.75; end; if carga==4 size_xlsx=size(P100,2);
84
P_IMEP=P100; BMEP=500*1.00; end;
for curva=2:size_xlsx
carga curva
if carga==1 P=P25(:,curva); end; if carga==2 P=P50(:,curva); end; if carga==3 P=P75(:,curva); end; if carga==4 P=P100(:,curva); end;
%% Selecionando os dados a trabalhar %DADOS AUXILIARES rc=15.2; %manual MAN B=160; L=240; %manual MAN a=L/2; %heywood p43 l=475.5; %medição da peça R=l/a; %heywood p43 vd=L*pi*(B/2)^2; vc=vd/(rc-1); %heywood p43
% abertura_valvula_admissao=-290; % tempo_de_abertura_valvula_admissao=320; % abertura_valvula_exaustao=90; inicio_compressao_isent=-75; final_compressao_isent=0; inicio_expansao_isent=0; final_expansao_isent=95;
angulo_a_mais_comp=50; angulo_a_mais_exp=50;
r2_minimo_comp=0.9999; r2_minimo_exp=0.9998; intervalo_angulo=0.25; %intervalo de angulo de coleta da curva de
pressão
pontos_min_para_ajuste_linear_comp=30/intervalo_angulo; %30 graus pontos_min_para_ajuste_linear_exp=30/intervalo_angulo; %30 graus ang_a_menos_para_segundo_ajuste_comp=10; ang_a_menos_para_segundo_ajuste_exp=10;
%% calculo e montagem da tabela do volume em funcao do angulo do
virabrequim ANGULO=-360-intervalo_angulo; for x=1:(720/intervalo_angulo+1) ANGULO=ANGULO+intervalo_angulo; %angulo de -360 a 360 com
intervalo de intervalo_angulo
85
%criacao da matriz do volume: coluna1=angulo; coluna2=volume da
camara %de combustao para cada ângulo v(x,1)=ANGULO; v(x,2)=vc*(1+.5*(rc-1)*(R+1-cos(degtorad(ANGULO))-(R^2-
(sin(degtorad(ANGULO)))^2)^.5)); %heywood p44
end;
%% Plotando o gráfico ângulo x volume % figure; % plot(v(:,1),v(:,2))
%% montagem da tabela da pressao em funcao do angulo do virabrequim ANGULO=-360-intervalo_angulo; for x=1:(720/intervalo_angulo+1) ANGULO=ANGULO+intervalo_angulo; pressao(x,1)=ANGULO; pressao(x,2)=P(x,1); end;
%% plotando o gráfico ângulo x pressão % figure; % plot(pressao(:,1),pressao(:,2))
%% montagem de duas matrizes: % 1) ln(volume) em funcao do angulo do virabrequim pra intervalo de % ângulo [-360,360] % 2) ln(pressao) em funcao do angulo do virabrequim pra intervalo de % ângulo [-360,360] ANGULO=-360-intervalo_angulo; for x=1:(720/intervalo_angulo+1) ANGULO=ANGULO+intervalo_angulo;
lnV(x,1)=ANGULO; lnV(x,2)=log(v(x,2)); %logarítimo neperiano
lnP(x,1)=ANGULO; lnP(x,2)=log(pressao(x,2));
end;
%% plotando gráfico lnV x lnP para ângulo de [-360,360] % figure; % plot(lnV(:,2),lnP(:,2)) % hold all
%% Compressao % ajuste do intervalo de apenas compressao lnV X lnP para uma funcao
linear tendo como critério uma reta com R²>r2_minimo_comp (valor
definido em Dados)
% Loop para definir melhor intervalo da compressão isentrópica sem
incluir pontos após o início da combustão
for ajuste_comp=1:2
if ajuste_comp==1
86
final_compressao_isent_1=final_compressao_isent+intervalo_angulo; %
variável nova que será, o valor do final da compressão isentrópica
obedecendo o critério do R² r2_comp=0; else % No loop em que ajuste_comp=2, o objetivo é diminuir alguns graus
e refazer outro ajuste linear. Se faz isso pois como há muitos pontos,
e o ponto do início da combustão não terá grande efeito no calculo de
R², permitindo que se consiga fazer o ajuste linear em um intervalo
contendo o início da combustão. Sendo assim terei um ajuste linear de
um intervalo de ângulo mesmo com um ponto bem fora do ajuste linear, o
que não é o desejável. length_lnV_comp=size(lnV_comp,1); if length_lnV_comp<pontos_min_para_ajuste_linear_comp %Se eu
tiver um número de pontos menor que um valor abitrário definido em
DADOS AUXILIARES (pontos_min_para_ajuste_linear_comp) eu não continuo
a calcular essa curva de pressão(apenas uma restrição para economizar
tempo já que a restrição se encontra após os cálculos).
r2_comp = r2_minimo_comp; else % Para refazer o calculo com alguns pontos a menos preciso
primeiro limpar algumas varáveis
varlist={'lnV_comp','lnP_comp','funcao_reta_comp','comp_reta_residuo'}
; clear(varlist{:}); final_compressao_isent_1=final_compressao_isent_1-
ang_a_menos_para_segundo_ajuste_comp+intervalo_angulo; r2_comp=0; end;
end;
%selecionando o intervalo para ajuste linear tendo como critério o R². while r2_comp < r2_minimo_comp
final_compressao_isent_1=final_compressao_isent_1-
intervalo_angulo; %diminui a cada iteração de intervalo_angulo ANGULO=inicio_compressao_isent-intervalo_angulo; %O intervalo diminui constantemente (mantem fixo o angulo de
início da compressão isentrópica, já definido, e aos poucos diminuo o
ângulo do final da compressão isentrópica) até encontrarmos o ajuste
linear para R² desejado
if abs(final_compressao_isent_1-inicio_compressao_isent)>0 % montando a matriz de lnV e lnP para um dado intervalo for x=1:abs(final_compressao_isent_1-
inicio_compressao_isent)/intervalo_angulo+1
ANGULO=ANGULO+intervalo_angulo; %converção do valor da
interação (x) para o ângulo que se está trabalhando. LINHA=1+(ANGULO+360)/intervalo_angulo; %Identificando a linha
na função lnV e lnP para o ângulo correspondente ao loop de x
lnV_comp(x,1)=ANGULO; lnV_comp(x,2)=lnV(LINHA,2); lnP_comp(x,1)=ANGULO;
87
lnP_comp(x,2)=lnP(LINHA,2); end;
%% Ajuste linear de uma reta na compressão para ln(Volume) e
ln(Pressão) funcao_reta_comp=fit(lnV_comp(:,2),lnP_comp(:,2),'poly1');
%% calculando o R^2 da aproximação linear
%Média de lnP soma_lnP=0; for x=1:abs(final_compressao_isent_1-
inicio_compressao_isent)/intervalo_angulo+1 soma_lnP=soma_lnP+lnP_comp(x,2); end; media_lnP=soma_lnP/x;
% Calculo da soma total do resíduo quadrado e soma total dos quadrados ANGULO=inicio_compressao_isent-intervalo_angulo; soma_total_dos_quadrados_comp=0; soma_do_quadrado_dos_residuos_comp=0;
for x=1:(final_compressao_isent_1-
inicio_compressao_isent)/intervalo_angulo+1
ANGULO=ANGULO+intervalo_angulo;
comp_reta_residuo(x,1)=ANGULO; comp_reta_residuo(x,2)=lnP_comp(x,2)-
funcao_reta_comp(lnV_comp(x,2));
soma_do_quadrado_dos_residuos_comp=soma_do_quadrado_dos_residuos_comp+
(comp_reta_residuo(x,2))^2;
soma_total_dos_quadrados_comp=soma_total_dos_quadrados_comp+(lnP_comp(
x,2)-media_lnP)^2;
end;
%cálculo de R² r2_comp=1-
soma_do_quadrado_dos_residuos_comp/soma_total_dos_quadrados_comp;
%limpando variáveis para refazer a conta se r2_minimo_comp não foi
atengido if r2_comp<r2_minimo_comp
varlist={'lnV_comp','lnP_comp','funcao_reta_comp','comp_reta_residuo'}
; clear(varlist{:}); end
else r2_comp=1.1; %parar cálculo end; end;
%mostrando r2_comp com 10 dígitos
88
vpa(r2_comp,10)
if r2_comp>1 lnV_comp=0; break % Para-se o "if" forçando o programa a ir para o "while",
que devido ao valor de r2_comp para o calculo end;
end;
% Incluindo critério que se tiver um número de pontos no novo ajuste
linear % menor que um valor arbritário não continuo com o calculo length_lnV_comp=size(lnV_comp,1);
if length_lnV_comp>pontos_min_para_ajuste_linear_comp %parar cálculo
dessa curva de pressão se esse critério for falso
% %% Plotando a reta na compressao dos valores de lnV para os angulos
de -180º a 0º % % ANGULO=-180-intervalo_angulo; % for x=1:(180/intervalo_angulo+1) % % ANGULO=ANGULO+intervalo_angulo; % LINHA=1+(ANGULO+360)/intervalo_angulo; % % %matriz com lnV; lnP ajuste linear; e ângulo correspondente % comp_reta_matriz(x,1)=lnV(LINHA,2); % comp_reta_matriz(x,2)=funcao_reta_comp(lnV(LINHA,2)); % comp_reta_matriz(x,3)=ANGULO; % % end; % % % Plotando ajuste linear na compressão % plot(comp_reta_matriz(:,1),comp_reta_matriz(:,2));
%% expansao % ajuste do intervalo de apenas expansao lnV X lnP para uma funcao
linear tendo como critério uma reta com R²>r2_minimo_exp (valor
definido em Dados)
% Loop para definir melhor intervalo da expressão isentrópica sem
incluir pontos após o início da combustão
for ajuste_exp=1:2
if ajuste_exp==1 inicio_expansao_isent_1=inicio_expansao_isent; % variável nova
que será, o valor do final da expressão isentrópica obedecendo o
critério do R² r2_exp=0; else % No loop em que ajuste_exp=2, o objetivo é diminuir alguns graus
e refazer outro ajuste linear. Se faz isso pois como há muitos pontos,
e o ponto do início da combustão não terá grande efeito no calculo de
R², permitindo que se consiga fazer o ajuste linear em um intervalo
contendo o início da combustão. Sendo assim terei um ajuste linear de
um intervalo de ângulo mesmo com um ponto bem fora do ajuste linear, o
que não é o desejável.
89
length_lnV_exp=size(lnV_exp,1); if length_lnV_exp<pontos_min_para_ajuste_linear_exp %Se eu
tiver um número de pontos menor que um valor abitrário definido em
DADOS AUXILIARES (pontos_min_para_ajuste_linear_exp) eu não continuo a
calcular essa curva de pressão(apenas uma restrição para economizar
tempo já que a restrição se encontra após os cálculos).
r2_exp=r2_minimo_exp; else % Para refazer o calculo com alguns pontos a menos preciso
primeiro limpar algumas varáveis
varlist={'lnV_exp','lnP_exp','funcao_reta_exp','exp_reta_residuo'}; clear(varlist{:});
inicio_expansao_isent_1=inicio_expansao_isent_1+ang_a_menos_para_segun
do_ajuste_exp+intervalo_angulo; r2_exp=0; end;
end;
%selecionando o intervalo para ajuste linear tendo como critério o R². while r2_exp < r2_minimo_exp
inicio_expansao_isent_1=inicio_expansao_isent_1+intervalo_angulo;
%diminui a cada iteração de intervalo_angulo ANGULO=inicio_expansao_isent_1-intervalo_angulo; %O intervalo diminui constantemente (mantem fixo o angulo de
início da expressão isentrópica, já definido, e aos poucos diminuo o
ângulo do final da expressão isentrópica) até encontrarmos o ajuste
linear para R² desejado
if abs(final_expansao_isent-inicio_expansao_isent_1)>0 % montando a matriz de lnV e lnP para um dado intervalo for x=1:abs(final_expansao_isent-
inicio_expansao_isent_1)/intervalo_angulo+1
ANGULO=ANGULO+intervalo_angulo; %converção do valor da
interação (x) para o ângulo que se está trabalhando. LINHA=1+(ANGULO+360)/intervalo_angulo; %Identificando a linha
na função lnV e lnP para o ângulo correspondente ao loop de x
lnV_exp(x,1)=ANGULO; lnV_exp(x,2)=lnV(LINHA,2); lnP_exp(x,1)=ANGULO; lnP_exp(x,2)=lnP(LINHA,2); end;
%% Ajuste linear de uma reta na expressão para ln(Volume) e
ln(Pressão) funcao_reta_exp=fit(lnV_exp(:,2),lnP_exp(:,2),'poly1');
%% calculando o R^2 da aproximação linear
%Média de lnP soma_lnP=0;
90
for x=1:abs(final_expansao_isent-
inicio_expansao_isent_1)/intervalo_angulo+1 soma_lnP=soma_lnP+lnP_exp(x,2); end; media_lnP=soma_lnP/x;
% Calculo da soma total do resíduo quadrado e soma total dos quadrados ANGULO=inicio_expansao_isent_1-intervalo_angulo; soma_total_dos_quadrados_exp=0; soma_do_quadrado_dos_residuos_exp=0;
for x=1:(final_expansao_isent-
inicio_expansao_isent_1)/intervalo_angulo+1
ANGULO=ANGULO+intervalo_angulo;
exp_reta_residuo(x,1)=ANGULO; exp_reta_residuo(x,2)=lnP_exp(x,2)-funcao_reta_exp(lnV_exp(x,2));
soma_do_quadrado_dos_residuos_exp=soma_do_quadrado_dos_residuos_exp+(e
xp_reta_residuo(x,2))^2;
soma_total_dos_quadrados_exp=soma_total_dos_quadrados_exp+(lnP_exp(x,2
)-media_lnP)^2;
end;
%cálculo de R² r2_exp=1-
soma_do_quadrado_dos_residuos_exp/soma_total_dos_quadrados_exp;
%limpando variáveis para refazer a conta se r2_minimo_exp não foi
atengido if r2_exp<r2_minimo_exp
varlist={'lnV_exp','lnP_exp','funcao_reta_exp','exp_reta_residuo'}; clear(varlist{:}); end
else r2_exp=1.1; %parar cálculo end; end;
%mostrando r2_exp com 10 dígitos vpa(r2_exp,10)
if r2_exp>1 lnV_exp=0; break % Para-se o "if" forçando o programa a ir para o "while",
que devido ao valor de r2_exp para o calculo end;
end;
% Incluindo critério que se tiver um número de pontos no novo ajuste
linear % menor que um valor arbritário não continuo com o calculo
91
length_lnV_exp=size(lnV_exp,1);
if length_lnV_exp>pontos_min_para_ajuste_linear_exp %parar cálculo
dessa curva de pressão se esse critério for falso
% %% Plotando a reta na expansao dos valores de lnV para os angulos de
-180º a 0º % % ANGULO=-180-intervalo_angulo; % for x=1:(180/intervalo_angulo+1) % % ANGULO=ANGULO+intervalo_angulo; % LINHA=1+(ANGULO+360)/intervalo_angulo; % % %matriz com lnV; lnP ajuste linear; e ângulo correspondente % exp_reta_matriz(x,1)=lnV(LINHA,2); % exp_reta_matriz(x,2)=funcao_reta_exp(lnV(LINHA,2)); % exp_reta_matriz(x,3)=ANGULO; % % end; % % % Plotando ajuste linear na expressão % plot(exp_reta_matriz(:,1),exp_reta_matriz(:,2));
%% Metodo das DERIVADAS do erro %% compressao
% matriz com a diferença de cada ponto entre o valor experimental de
lnP e o valor de lnP ajustado pela reta
ANGULO=inicio_compressao_isent-intervalo_angulo; for x=1:abs(inicio_compressao_isent-
(final_compressao_isent_1+angulo_a_mais_comp))/intervalo_angulo
ANGULO=ANGULO+intervalo_angulo; LINHA=1+(ANGULO+360)/intervalo_angulo;
erro_comp(x,1)=ANGULO; erro_comp(x,2)=(lnP(LINHA,2)-funcao_reta_comp(lnV(LINHA,2)));
end;
% % Plotando % figure; % plot(erro_comp(:,1),erro_comp(:,2))
%% Calculando a derivada do erro_comp. São criadas duas matrizes: 1ª
número de pontos igual ao erro_comp e a 2ª apenas com os pontos
responsáveis por gerar o ajuste linear com critério definido de um R²
mínimo ANGULO=inicio_compressao_isent; for g=1:2;
if g==1 H=final_compressao_isent_1+angulo_a_mais_comp-
intervalo_angulo; else H=final_compressao_isent_1-intervalo_angulo;
92
end;
for x=1:abs(inicio_compressao_isent-H)/intervalo_angulo
if g==1 % Calculo aproximado da derivada de erro_comp. Sendo o calculo
feito entre cada dois pontos. O calculo é feito para o mesmo intervalo
de erro_comp. der_comp(x,1)=(erro_comp(x+1,1)-
erro_comp(x,1))/2+erro_comp(x,1); der_comp(x,2)=(erro_comp(x+1,2)-
erro_comp(x,2))/(erro_comp(x+1,1)-erro_comp(x,1)); end;
if g==2 %Calculo aproximado da derivada do erro_comp entre cada dois
pontos para o intervalo de angulos em que se faz o ajuste linear para
a compressão (R²>r2_minimo_comp) der_comp_1(x,1)=(erro_comp(x+1,1)-
erro_comp(x,1))/2+erro_comp(x,1); der_comp_1(x,2)=(erro_comp(x+1,2)-
erro_comp(x,2))/(erro_comp(x+1,1)-erro_comp(x,1)); g2=x; end; end; end;
%% encontrando o maximo valor absoluto da der_comp_1, ou seja, para a
derivada do intervalo de angulo que se conseguiu fazer um ajuste
linear com R²>r2_minimo_comp max_abs_der_comp=abs(der_comp_1(1,2)); for y=2:g2 if abs(der_comp_1(y,2))>max_abs_der_comp max_abs_der_comp=abs(der_comp_1(y,2)); end; end;
% % Plotando a der_comp para um intervalo de ângulo maior que o do
ajuste linear da compressão (R²>r2_minimo_comp) % figure; % plot(der_comp(:,1),der_comp(:,2)) % hold all % %A partir do valor máximo absoluto da função dif_der_comp, plota-se
duas % %retas paralelas ao eixo y, uma com esse valor positivo e outra com
o valor % %negativo % plot(der_comp(:,1),max_abs_der_comp); % plot(der_comp(:,1),-max_abs_der_comp); % hold off
%% encontrando o angulo de inicio da combustão
%diferença de cada ponto entre o valor experimental de lnP e o valor
de lnP ajustado pela reta
93
% O critério utilizado foi o ângulo em que a der_comp se tornou 10%
maior que o maior der_comp encontrada no intervalo % em que se conseguiu fazer o ajuste linear da compressão
(R²>r2_minimo_comp) y=2; b=0; length_comp=size(der_comp,1); while abs(b)<1.5*max_abs_der_comp y=y+1; if y<=length_comp-2 ang=der_comp(y,1); b=der_comp(y,2); % Como a turbulência gerada pela injeção de combustível afeta o
valor de erro_comp, fazendo este ultimo oscilar bastante, então, para
não encontrar erroneamente o início da combustão será posta mais uma
restrição. A restrição é que a média de alguns pontos de erro_comp
deve ser maior que o valor do maior erro_comp no intervalo usado para
ajustar a reta. Como o salto dado pelo início da combustão é muito
grande e os seguintes também, a restrição a seguir será usada if abs(b)>1.5*max_abs_der_comp if (der_comp(y-
1,2)+der_comp(y,2)+der_comp(y+1,2)+der_comp(y+2,2))/4<2.5*max_abs_der_
comp b=max_abs_der_comp/2; end; end;
else b=999*max_abs_der_comp; end; end;
if y<=length_comp-2 angulo_inicio_da_combustao=ang-intervalo_angulo/2
%% Metodo das DERIVADAS do erro %% expansao
% matriz com a diferença de cada ponto entre o valor experimental de
lnP e o valor de lnP ajustado pela reta
ANGULO=(inicio_expansao_isent_1-angulo_a_mais_exp); for x=1:abs((inicio_expansao_isent_1-angulo_a_mais_exp)-
final_expansao_isent)/intervalo_angulo
ANGULO=ANGULO+intervalo_angulo; LINHA=1+(ANGULO+360)/intervalo_angulo;
erro_exp(x,1)=ANGULO; erro_exp(x,2)=(lnP(LINHA,2)-funcao_reta_exp(lnV(LINHA,2)));
end;
% % Plotando % figure; % plot(erro_exp(:,1),erro_exp(:,2))
%% Calculando a derivada do erro_exp. São criadas duas matrizes: 1ª
número de pontos igual ao erro_exp e a 2ª apenas com os pontos
94
responsáveis por gerar o ajuste linear com critério definido de um R²
mínimo ANGULO=inicio_expansao_isent_1; for g=1:2;
if g==1 H=inicio_expansao_isent_1-angulo_a_mais_exp; else H=inicio_expansao_isent_1; end;
for x1=1:abs(H-final_expansao_isent)/intervalo_angulo-1 x=size(erro_exp,1)-x1;
if g==1 % Calculo aproximado da derivada de erro_exp. Sendo o calculo
feito entre cada dois pontos. O calculo é feito para o mesmo intervalo
de erro_exp. % der_exp(x,1)=(erro_exp(x+1,1)-
erro_exp(x,1))/2+erro_exp(x,1); % der_exp(x,2)=(erro_exp(x+1,2)-
erro_exp(x,2))/(erro_exp(x+1,1)-erro_exp(x,1)); erro_exp_2(x,1)=erro_exp(x,1); erro_exp_2(x,2)=erro_exp(x,2); end;
if g==2 %Calculo aproximado da derivada do erro_exp entre cada dois
pontos para o intervalo de angulos em que se faz o ajuste linear para
a expansão (R²>r2_minimo_exp) % der_exp_1(x,1)=(erro_exp(x+1,1)-
erro_exp(x,1))/2+erro_exp(x,1); % der_exp_1(x,2)=(erro_exp(x+1,2)-
erro_exp(x,2))/(erro_exp(x+1,1)-erro_exp(x,1)); erro_exp_2_1(x,1)=erro_exp(x,1); erro_exp_2_1(x,2)=erro_exp(x,2);
end; end; end; x_exp_2_1=x;
%% encontrando o maximo valor absoluto da der_exp_1, ou seja, para a
derivada do intervalo de angulo que se conseguiu fazer um ajuste
linear com R²>r2_minimo_exp max_abs_erro_2_exp=abs(erro_exp_2_1(1,2)); for y=2:size(erro_exp_2_1,1) if abs(erro_exp_2_1(y,2))>max_abs_erro_2_exp max_abs_erro_2_exp=abs(erro_exp_2_1(y,2)); end; end;
% % Plotando a der_exp para um intervalo de ângulo maior que o do
ajuste linear da expansão (R²>r2_minimo_exp) % figure; % plot(der_exp(:,1),der_exp(:,2)); % hold all % %A partir do valor máximo absoluto da função dif_der_exp, plota-se
duas
95
% %retas paralelas ao eixo y, uma com esse valor positivo e outra com
o valor % %negativo % plot(der_exp(:,1),max_abs_der_exp); % plot(der_exp(:,1),-max_abs_der_exp); % hold off
%% encontrando o angulo de inicio da combustão
%diferença de cada ponto entre o valor experimental de lnP e o valor
de lnP ajustado pela reta
% O critério utilizado foi o ângulo em que a der_exp se tornou 10%
maior que o maior der_exp encontrada no intervalo % em que se conseguiu fazer o ajuste linear da expansão
(R²>r2_minimo_exp)
b=0; length_exp=size(erro_exp_2,1); y=length_exp-2; while abs(b)<1.02*max_abs_erro_2_exp y=y-1; if y>3 ang=erro_exp_2(y,1); b=erro_exp_2(y,2);
% Como a turbulência gerada pela injeção de combustível afeta o
valor de erro_exp, fazendo este ultimo oscilar bastante, então, para
não encontrar erroneamente o início da combustão será posta mais uma
restrição. A restrição é que a média de alguns pontos de erro_exp deve
ser maior que o valor do maior erro_exp no intervalo usado para
ajustar a reta. Como o salto dado pelo início da combustão é muito
grande e os seguintes também, a restrição a seguir será usada % garantir que todos os valores antes do ponto que se esta
analisando estão abaixo da linha da regressão linear if y<length_exp-85 if y>80 n_der_erro_menor=0; for d=1:80 if abs(erro_exp_2(y-
d,2))<1.01*max_abs_erro_2_exp; n_der_erro_menor=n_der_erro_menor+1; end; end; sempre_negativos_abs=0; sempre_negativos=0; for c=1:80 sempre_negativos_abs=abs(erro_exp(y-
c,2))+sempre_negativos_abs; sempre_negativos=erro_exp(y-
c,2)+sempre_negativos; end; sempre_negativos=abs(sempre_negativos); if sempre_negativos~=sempre_negativos_abs b=0.1*max_abs_erro_2_exp; end; if n_der_erro_menor>0 b=0.1*max_abs_erro_2_exp; end; end;
96
end;
else b=999*max_abs_der_exp; end; end; if y>pontos_min_para_ajuste_linear_exp angulo_final_da_combustao=ang-intervalo_angulo/2 if y>x_exp_2_1;
%% IMEP
%% IMEP _ Trapezoidal rule
for y=1:4 IMEP_matriz(y)=0; for x=(y-1)*180/intervalo_angulo+1:y*180/intervalo_angulo
IMEP_matriz(y)=IMEP_matriz(y)+0.5*(P(x,1)+P(x+1,1))*abs(v(x,2)-
v(x+1,2)); end; end;
IMEP_net=(IMEP_matriz(3)-IMEP_matriz(2))-(IMEP_matriz(4)-
IMEP_matriz(1));
BMEP;
FMEP=IMEP_net-BMEP; % Heywod página 47 e 48
%% Pmax e seu respectivo ângulo
ANGULO=-360-intervalo_angulo; Pmax=P(1); for x=2:(720/intervalo_angulo+1) ANGULO=ANGULO+intervalo_angulo; %angulo de -360 a 360 com
intervalo de intervalo_angulo if P(x)>Pmax Pmax=P(x); ANGmax=ANGULO; end; end; Pmax ANGmax
%% montar uma matriz com os valores calculados do inicio e do final da
combustão para cada curva de pressão que se conseguiu calcular esse
valor pelo método utilizado q=q+1; matriz(q,1)=curva; matriz(q,2)=angulo_inicio_da_combustao;
%degree matriz(q,3)=angulo_final_da_combustao;
%degree matriz(q,4)=IMEP_net*5/10^6;
%kW matriz(q,5)=BMEP;
%kW
97
matriz(q,6)=matriz(q,4)-matriz(q,5);
%kW FMEP matriz(q,7)=BMEP/(((10^5*10^-9*(1200/2)*(1/60)*5)/10^3)*vd);
%bar matriz(q,8)=matriz(q,4)/(((10^5*10^-9*(1200/2)*(1/60)*5)/10^3)*vd);
%bar matriz(q,9)=matriz(q,8)-matriz(q,7);
%bar FMEP matriz(q,10)=Pmax; matriz(q,11)=ANGmax;
%% end; end; end; end; end;
clearvars -except BMEP curva matriz P q P100 P75 P50 P25 matriz_P25
matriz_P50 matriz_P75 matriz_P100 carga size_xlsx dados_validos_P25
dados_validos_P50 dados_validos_P75 dados_validos_P100 P_IMEP
intervalo_angulo media_desvPad_matriz_P25 media_desvPad_matriz_P50
media_desvPad_matriz_P75 media_desvPad_matriz_P100
end;
%% Exportando os dados da combustão para um arquido do Excel if q>0; if carga==1 matriz_P25=matriz; xlswrite('matriz angulos cada curva',matriz_P25,'P25'); %Exportando
os dados encontrados para um arquivo do excel end;
if carga==2 matriz_P50=matriz; xlswrite('matriz angulos cada curva',matriz_P50,'P50'); %Exportando
os dados encontrados para um arquivo do excel end;
if carga==3 matriz_P75=matriz; xlswrite('matriz angulos cada curva',matriz_P75,'P75'); %Exportando
os dados encontrados para um arquivo do excel end;
if carga==4 matriz_P100=matriz; xlswrite('matriz angulos cada curva',matriz_P100,'P100'); %Exportando
os dados encontrados para um arquivo do excel end; end;
clear matriz;
98
%% FASE DOIS --> ESCOLHENDO OS DADOS "BONS"
%% Escolhendo a matriz a se trabalhar
varlist={'matriz_P','n','media_matriz_P', 'standard_deviation',
'chauvenet_criterion', 'chauvenet_solve', 'chauvenet_inicio_comb',
'chauvenet_final_comb'}; clear(varlist{:});
if carga==1 matriz_P=matriz_P25; end; if carga==2 matriz_P=matriz_P50; end; if carga==3 matriz_P=matriz_P75; end; if carga==4 matriz_P=matriz_P100; end;
%% Calculando a média do incicio da combustão (coluna 1), média do
final da combustão (coluna 2), média do final da combustão pelo
metodo2 (coluna 3)
n=size(matriz_P,1); for y=1:size(matriz_P,2)-1 soma_matriz_P(y)=0; for x=1:n soma_matriz_P(y)=soma_matriz_P(y)+matriz_P(x,y+1); end; media_matriz_P(y)=soma_matriz_P(y)/x;
%% % Calculando a o Standard Deviation standard_deviation(y)=0; for x=1:n standard_deviation(y)=standard_deviation(y)+(matriz_P(x,y+1)-
media_matriz_P(1,y))^2; end; standard_deviation(y)=sqrt(standard_deviation(y)/(n-1));
end;
%% Calculando o valor máximo para (deviation/Satandard_Deviation)
baseando-se no critério de Chauvenet. Critério de Chauvenet pode ser
encontrado no livro: Holman, J. P.; Experimental Methods for
engineers, 6th edition; Página 70) syms chauvenet
chauvenet_solve=solve(0.5*erf(chauvenet/sqrt(2))==(1-
1/(2*n))/2,chauvenet); chauvenet_criterion=vpa(chauvenet_solve);
99
%% criar uma matriz apenas com os valores que foram considerados
validos simultaneamente no final e no inicio da combustão z=0; for x=1:n chauvenet_inicio_comb=abs((matriz_P(x,2)-
media_matriz_P(1))/standard_deviation(1)); chauvenet_final_comb=abs((matriz_P(x,3)-
media_matriz_P(2))/standard_deviation(2)); chauvenet_IMEP=abs((matriz_P(x,4)-
media_matriz_P(3))/standard_deviation(3));
if chauvenet_inicio_comb<chauvenet_criterion if chauvenet_final_comb<chauvenet_criterion if chauvenet_IMEP<chauvenet_criterion z=z+1; if carga==1 for x1=1:size(matriz_P,2) dados_validos_P25(z,x1)=matriz_P(x,x1); end; end; if carga==2 for x1=1:size(matriz_P,2) dados_validos_P50(z,x1)=matriz_P(x,x1); end; end; if carga==3 for x1=1:size(matriz_P,2) dados_validos_P75(z,x1)=matriz_P(x,x1); end; end; if carga==4 for x1=1:size(matriz_P,2) dados_validos_P100(z,x1)=matriz_P(x,x1); end; end; end; end; end; end;
if z>0; if carga==1 n=size(dados_validos_P25,1); for y=1:size(dados_validos_P25,2)-1 soma_matriz_P25(y)=0; for x=1:n
soma_matriz_P25(y)=soma_matriz_P25(y)+dados_validos_P25(x,y+1); end; media_desvPad_matriz_P25(y,1)=soma_matriz_P25(y)/x;
media_desvPad_matriz_P25(y,2)=0; for x=1:n
media_desvPad_matriz_P25(y,2)=media_desvPad_matriz_P25(y,2)+(dados_val
idos_P25(x,y+1)-media_desvPad_matriz_P25(y,1))^2; end;
media_desvPad_matriz_P25(y,2)=sqrt(media_desvPad_matriz_P25(y,2)/(n-
1)); end;
100
media_desvPad_matriz_P25(1,3)=size(matriz_P25,1); media_desvPad_matriz_P25(2,3)=size(dados_validos_P25,1); xlswrite('Resultado.xls',media_desvPad_matriz_P25,'P25'); end;
if carga==2 n=size(dados_validos_P50,1); for y=1:size(dados_validos_P50,2)-1 soma_matriz_P50(y)=0; for x=1:n
soma_matriz_P50(y)=soma_matriz_P50(y)+dados_validos_P50(x,y+1); end; media_desvPad_matriz_P50(y,1)=soma_matriz_P50(y)/x;
media_desvPad_matriz_P50(y,2)=0; for x=1:n
media_desvPad_matriz_P50(y,2)=media_desvPad_matriz_P50(y,2)+(dados_val
idos_P50(x,y+1)-media_desvPad_matriz_P50(y,1))^2; end;
media_desvPad_matriz_P50(y,2)=sqrt(media_desvPad_matriz_P50(y,2)/(n-
1)); end; media_desvPad_matriz_P50(1,3)=size(matriz_P50,1); media_desvPad_matriz_P50(2,3)=size(dados_validos_P50,1); xlswrite('Resultado.xls',media_desvPad_matriz_P50,'P50'); end;
if carga==3 n=size(dados_validos_P75,1); for y=1:size(dados_validos_P75,2)-1 soma_matriz_P75(y)=0; for x=1:n
soma_matriz_P75(y)=soma_matriz_P75(y)+dados_validos_P75(x,y+1); end; media_desvPad_matriz_P75(y,1)=soma_matriz_P75(y)/x;
media_desvPad_matriz_P75(y,2)=0; for x=1:n
media_desvPad_matriz_P75(y,2)=media_desvPad_matriz_P75(y,2)+(dados_val
idos_P75(x,y+1)-media_desvPad_matriz_P75(y,1))^2; end;
media_desvPad_matriz_P75(y,2)=sqrt(media_desvPad_matriz_P75(y,2)/(n-
1)); end; media_desvPad_matriz_P75(1,3)=size(matriz_P75,1); media_desvPad_matriz_P75(2,3)=size(dados_validos_P75,1); xlswrite('Resultado.xls',media_desvPad_matriz_P75,'P75'); end;
if carga==4 n=size(dados_validos_P100,1); for y=1:size(dados_validos_P100,2)-1 soma_matriz_P100(y)=0; for x=1:n
101
soma_matriz_P100(y)=soma_matriz_P100(y)+dados_validos_P100(x,y+1); end; media_desvPad_matriz_P100(y,1)=soma_matriz_P100(y)/x;
media_desvPad_matriz_P100(y,2)=0; for x=1:n
media_desvPad_matriz_P100(y,2)=media_desvPad_matriz_P100(y,2)+(dados_v
alidos_P100(x,y+1)-media_desvPad_matriz_P100(y,1))^2; end;
media_desvPad_matriz_P100(y,2)=sqrt(media_desvPad_matriz_P100(y,2)/(n-
1)); end; media_desvPad_matriz_P100(1,3)=size(matriz_P100,1); media_desvPad_matriz_P100(2,3)=size(dados_validos_P100,1); xlswrite('Resultado.xls',media_desvPad_matriz_P100,'P100'); end; end;
varlist={'soma_matriz_P25','soma_matriz_P50','soma_matriz_P75','soma_m
atriz_P100'}; clear(varlist{:});
end;
%% Encontrando curva representativa
for carga=1:4
clear dados
if carga==1 P=P25(:,curva); dados=dados_validos_P25; dados_final=media_desvPad_matriz_P25; end; if carga==2 P=P50(:,curva); dados=dados_validos_P50; dados_final=media_desvPad_matriz_P50; end; if carga==3 P=P75(:,curva); dados=dados_validos_P75; dados_final=media_desvPad_matriz_P75; end; if carga==4 P=P100(:,curva); dados=dados_validos_P100; dados_final=media_desvPad_matriz_P100; end;
clear comparacao
for x=1:size(dados,1) comparacao(x,1)=dados(x,1); comparacao(x,2)=(dados(x,10)-dados_final(9,1))/dados(x,10); comparacao(x,3)=(dados(x,11)-dados_final(10,1))/dados(x,11);
102
comparacao(x,4)=(abs(2*comparacao(x,2))+abs(comparacao(x,3)))/3; end;
aux_comparacao(1)=comparacao(1,1); aux_comparacao(2)=comparacao(1,4); for x=2:size(comparacao,1) if comparacao(x,4)<aux_comparacao(2); aux_comparacao(1)=comparacao(x,1); aux_comparacao(2)=comparacao(x,4); end; end;
aux_comparacao
curva_padrao_1(:,1)=P25(:,1); %Angulos
for x=1:4 if carga==x curva_padrao(x)=aux_comparacao(1); curva_padrao_1(:,(x+1))=P(:,1); end; end;
end;
xlswrite('curvas.xls',curva_padrao_1);
103
Anexo II: Programa para calcular dados da curva de pressão e utilizando o método
descrito no item 4.5.2.2.
%% Limpando todas as variaveis e importando os dados clear;
% Importando as curvas de pressão de três tomadas diferentes e em
seguida juntando as curvas de mesma carga mas de tomadas diferentes em
uma única matriz
P25_1=importdata('Tomada 01\P25.txt'); P25_1=P25_1.data; P25_2=importdata('Tomada 02\P25.txt'); P25_2=P25_2.data; P25_3=importdata('Tomada 03\P25.txt'); P25_3=P25_3.data;
P50_1=importdata('Tomada 01\P50.txt'); P50_1=P50_1.data; P50_2=importdata('Tomada 02\P50.txt'); P50_2=P50_2.data; P50_3=importdata('Tomada 03\P50.txt'); P50_3=P50_3.data;
P75_1=importdata('Tomada 01\P75.txt'); P75_1=P75_1.data; P75_2=importdata('Tomada 02\P75.txt'); P75_2=P75_2.data; P75_3=importdata('Tomada 03\P75.txt'); P75_3=P75_3.data;
P100_1=importdata('Tomada 01\P100.txt'); P100_1=P100_1.data; P100_2=importdata('Tomada 02\P100.txt'); P100_2=P100_2.data; P100_3=importdata('Tomada 03\P100.txt'); P100_3=P100_3.data;
for y=1:3 if y==1; P25=P25_1; else if y==2 for x=2:size(P25_2,2) P25(:,size(P25,2)+1)=P25_2(:,x); end; else for x=2:size(P25_3,2) P25(:,size(P25,2)+1)=P25_3(:,x); end; end end
if y==1; P50=P50_1; else if y==2 for x=2:size(P50_2,2)
104
P50(:,size(P50,2)+1)=P50_2(:,x); end; else for x=2:size(P50_3,2) P50(:,size(P50,2)+1)=P50_3(:,x); end; end end
if y==1; P75=P75_1; else if y==2 for x=2:size(P75_2,2) P75(:,size(P75,2)+1)=P75_2(:,x); end; else for x=2:size(P75_3,2) P75(:,size(P75,2)+1)=P75_3(:,x); end; end end
if y==1; P100=P100_1; else if y==2 for x=2:size(P100_2,2) P100(:,size(P100,2)+1)=P100_2(:,x); end; else for x=2:size(P100_3,2) P100(:,size(P100,2)+1)=P100_3(:,x); end; end end
end;
%% Loop para calcular as curvas de pressão de diferentes cargas for carga=1:4 q=0;
if carga==1 size_xlsx=size(P25,2); P_IMEP=P25; BMEP=500*0.25; end; if carga==2 size_xlsx=size(P50,2); P_IMEP=P50; BMEP=500*0.50; end; if carga==3 size_xlsx=size(P75,2); P_IMEP=P75; BMEP=500*0.75; end; if carga==4 size_xlsx=size(P100,2);
105
P_IMEP=P100; BMEP=500*1.00; end;
for curva=2:size_xlsx
carga curva
if carga==1 P=P25(:,curva); end; if carga==2 P=P50(:,curva); end; if carga==3 P=P75(:,curva); end; if carga==4 P=P100(:,curva); end;
%% Selecionando os dados a trabalhar %DADOS AUXILIARES rc=15.2; %manual MAN B=160; L=240; %manual MAN a=L/2; %heywood p43 l=475.5; %medição da peça R=l/a; %heywood p43 vd=L*pi*(B/2)^2; vc=vd/(rc-1); %heywood p43
% abertura_valvula_admissao=-290; % tempo_de_abertura_valvula_admissao=320; % abertura_valvula_exaustao=90; inicio_compressao_isent=-75; final_compressao_isent=0; inicio_expansao_isent=0; final_expansao_isent=95;
angulo_a_mais_comp=50; angulo_a_mais_exp=120;
r2_minimo_comp=0.9999; r2_minimo_exp=0.9998; intervalo_angulo=0.25; %intervalo de angulo de coleta da curva de
pressão
pontos_min_para_ajuste_linear_comp=30/intervalo_angulo; %30 graus pontos_min_para_ajuste_linear_exp=30/intervalo_angulo; %30 graus ang_a_menos_para_segundo_ajuste_comp=10;
%% calculo e montagem da tabela do volume em funcao do angulo do
virabrequim ANGULO=-360-intervalo_angulo; for x=1:(720/intervalo_angulo+1)
106
ANGULO=ANGULO+intervalo_angulo; %angulo de -360 a 360 com
intervalo de intervalo_angulo %criacao da matriz do volume: coluna1=angulo; coluna2=volume da
camara %de combustao para cada ângulo v(x,1)=ANGULO; v(x,2)=vc*(1+.5*(rc-1)*(R+1-cos(degtorad(ANGULO))-(R^2-
(sin(degtorad(ANGULO)))^2)^.5)); %heywood p44
end;
%% Plotando o gráfico ângulo x volume % figure; % plot(v(:,1),v(:,2))
%% montagem da tabela da pressao em funcao do angulo do virabrequim ANGULO=-360-intervalo_angulo; for x=1:(720/intervalo_angulo+1) ANGULO=ANGULO+intervalo_angulo; pressao(x,1)=ANGULO; pressao(x,2)=P(x,1); end;
%% plotando o gráfico ângulo x pressão % figure; % plot(pressao(:,1),pressao(:,2))
%% montagem de duas matrizes: % 1) ln(volume) em funcao do angulo do virabrequim pra intervalo de % ângulo [-360,360] % 2) ln(pressao) em funcao do angulo do virabrequim pra intervalo de % ângulo [-360,360] ANGULO=-360-intervalo_angulo; for x=1:(720/intervalo_angulo+1) ANGULO=ANGULO+intervalo_angulo;
lnV(x,1)=ANGULO; lnV(x,2)=log(v(x,2)); %logarítimo neperiano
lnP(x,1)=ANGULO; lnP(x,2)=log(pressao(x,2));
end;
%% plotando gráfico lnV x lnP para ângulo de [-360,360] % figure; % plot(lnV(:,2),lnP(:,2)) % hold all
%% Compressao % ajuste do intervalo de apenas compressao lnV X lnP para uma funcao
linear tendo como critério uma reta com R²>r2_minimo_comp (valor
definido em Dados)
% Loop para definir melhor intervalo da compressão isentrópica sem
incluir pontos após o início da combustão
for ajuste_comp=1:2
107
if ajuste_comp==1
final_compressao_isent_1=final_compressao_isent+intervalo_angulo; %
variável nova que será, o valor do final da compressão isentrópica
obedecendo o critério do R² r2_comp=0; else % No loop em que ajuste_comp=2, o objetivo é diminuir alguns graus
e refazer outro ajuste linear. Se faz isso pois como há muitos pontos,
e o ponto do início da combustão não terá grande efeito no calculo de
R², permitindo que se consiga fazer o ajuste linear em um intervalo
contendo o início da combustão. Sendo assim terei um ajuste linear de
um intervalo de ângulo mesmo com um ponto bem fora do ajuste linear, o
que não é o desejável. length_lnV_comp=size(lnV_comp,1); if length_lnV_comp<pontos_min_para_ajuste_linear_comp %Se eu
tiver um número de pontos menor que um valor abitrário definido em
DADOS AUXILIARES (pontos_min_para_ajuste_linear_comp) eu não continuo
a calcular essa curva de pressão(apenas uma restrição para economizar
tempo já que a restrição se encontra após os cálculos).
r2_comp = r2_minimo_comp; else % Para refazer o calculo com alguns pontos a menos preciso
primeiro limpar algumas varáveis
varlist={'lnV_comp','lnP_comp','funcao_reta_comp','comp_reta_residuo'}
; clear(varlist{:}); final_compressao_isent_1=final_compressao_isent_1-
ang_a_menos_para_segundo_ajuste_comp+intervalo_angulo; r2_comp=0; end;
end;
%selecionando o intervalo para ajuste linear tendo como critério o R². while r2_comp < r2_minimo_comp
final_compressao_isent_1=final_compressao_isent_1-
intervalo_angulo; %diminui a cada iteração de intervalo_angulo ANGULO=inicio_compressao_isent-intervalo_angulo; %O intervalo diminui constantemente (mantem fixo o angulo de
início da compressão isentrópica, já definido, e aos poucos diminuo o
ângulo do final da compressão isentrópica) até encontrarmos o ajuste
linear para R² desejado
if abs(final_compressao_isent_1-inicio_compressao_isent)>0 % montando a matriz de lnV e lnP para um dado intervalo for x=1:abs(final_compressao_isent_1-
inicio_compressao_isent)/intervalo_angulo+1
ANGULO=ANGULO+intervalo_angulo; %converção do valor da
interação (x) para o ângulo que se está trabalhando. LINHA=1+(ANGULO+360)/intervalo_angulo; %Identificando a linha
na função lnV e lnP para o ângulo correspondente ao loop de x
lnV_comp(x,1)=ANGULO; lnV_comp(x,2)=lnV(LINHA,2);
108
lnP_comp(x,1)=ANGULO; lnP_comp(x,2)=lnP(LINHA,2); end;
%% Ajuste linear de uma reta na compressão para ln(Volume) e
ln(Pressão) funcao_reta_comp=fit(lnV_comp(:,2),lnP_comp(:,2),'poly1');
%% calculando o R^2 da aproximação linear
%Média de lnP soma_lnP=0; for x=1:abs(final_compressao_isent_1-
inicio_compressao_isent)/intervalo_angulo+1 soma_lnP=soma_lnP+lnP_comp(x,2); end; media_lnP=soma_lnP/x;
% Calculo da soma total do resíduo quadrado e soma total dos quadrados ANGULO=inicio_compressao_isent-intervalo_angulo; soma_total_dos_quadrados_comp=0; soma_do_quadrado_dos_residuos_comp=0;
for x=1:(final_compressao_isent_1-
inicio_compressao_isent)/intervalo_angulo+1
ANGULO=ANGULO+intervalo_angulo;
comp_reta_residuo(x,1)=ANGULO; comp_reta_residuo(x,2)=lnP_comp(x,2)-
funcao_reta_comp(lnV_comp(x,2));
soma_do_quadrado_dos_residuos_comp=soma_do_quadrado_dos_residuos_comp+
(comp_reta_residuo(x,2))^2;
soma_total_dos_quadrados_comp=soma_total_dos_quadrados_comp+(lnP_comp(
x,2)-media_lnP)^2;
end;
%cálculo de R² r2_comp=1-
soma_do_quadrado_dos_residuos_comp/soma_total_dos_quadrados_comp;
%limpando variáveis para refazer a conta se r2_minimo_comp não foi
atengido if r2_comp<r2_minimo_comp
varlist={'lnV_comp','lnP_comp','funcao_reta_comp','comp_reta_residuo'}
; clear(varlist{:}); end
else r2_comp=1.1; %parar cálculo end; end;
109
%mostrando r2_comp com 10 dígitos vpa(r2_comp,10)
if r2_comp>1 lnV_comp=0; break % Para-se o "if" forçando o programa a ir para o "while",
que devido ao valor de r2_comp para o calculo end;
end;
% Incluindo critério que se tiver um número de pontos no novo ajuste
linear % menor que um valor arbritário não continuo com o calculo length_lnV_comp=size(lnV_comp,1);
if length_lnV_comp>pontos_min_para_ajuste_linear_comp %parar cálculo
dessa curva de pressão se esse critério for falso
% %% Plotando a reta na compressao dos valores de lnV para os angulos
de -180º a 0º % % ANGULO=-180-intervalo_angulo; % for x=1:(180/intervalo_angulo+1) % % ANGULO=ANGULO+intervalo_angulo; % LINHA=1+(ANGULO+360)/intervalo_angulo; % % %matriz com lnV; lnP ajuste linear; e ângulo correspondente % comp_reta_matriz(x,1)=lnV(LINHA,2); % comp_reta_matriz(x,2)=funcao_reta_comp(lnV(LINHA,2)); % comp_reta_matriz(x,3)=ANGULO; % % end; % % Plotando ajuste linear na compressão % plot(comp_reta_matriz(:,1),comp_reta_matriz(:,2));
%% método para encontrar o final da combustão -> menor valor de
d(lnP)/d(teta) (derivada de lnP em função do ângulo)
ANGULO=0-intervalo_angulo; x_FC=1; der_exp(x_FC,2)=0; for x=1:abs(ANGULO-final_expansao_isent)/intervalo_angulo ANGULO=ANGULO+intervalo_angulo; LINHA=1+(ANGULO+360)/intervalo_angulo; % Calculo aproximado da derivada de erro_comp. Sendo ocalculo
feito entre cada dois pontos. der_exp(x,1)=(lnP(LINHA+1,1)-lnP(LINHA,1))/2+lnP(LINHA,1); der_exp(x,2)=(lnP(LINHA+1,2)-lnP(LINHA,2))/(lnP(LINHA+1,1)-
lnP(LINHA,1)); end;
ANGULO=0-intervalo_angulo;
menores_derivadas_0=(der_exp(1,2)+der_exp(2,2)+der_exp(3,2)+der_exp(4,
2)+der_exp(5,2)+der_exp(6,2)+der_exp(7,2)); for x=4:abs(ANGULO-final_expansao_isent)/intervalo_angulo-4
110
menores_derivadas=(der_exp(x-3,2)+der_exp(x-2,2)+der_exp(x-
1,2)+der_exp(x,2)+der_exp(x+1,2)+der_exp(x+2,2)+der_exp(x+3,2)); if menores_derivadas<menores_derivadas_0 menores_derivadas_0=menores_derivadas; x_FC=x; end; end;
intervalo=abs(7.5-der_exp(x_FC,1)-(intervalo_angulo/2)); if (intervalo+der_exp(x_FC,1)-
(intervalo_angulo/2))>final_expansao_isent intervalo=final_expansao_isent-der_exp(x_FC,1)-
(intervalo_angulo/2); end;
%Criando uma matriz com um intervalo de ângulo em torno do ponto %identificado por menores_derivadas aux1=0; for x=x_FC-
intervalo/intervalo_angulo:x_FC+intervalo/intervalo_angulo aux1=aux1+1; der_exp_min(aux1,1)=der_exp(x,1); der_exp_min(aux1,2)=der_exp(x,2); end;
% Aproximando os pontos para uma função polinomiala fim de
minimizar a influência do ruido dos dados experimentais funcao_der_exp_min=fit(der_exp_min(:,1),der_exp_min(:,2),'poly3');
% encontrando o ângulo em que a função polinomial tem seu mínimo
dentro do intervalo usado para cria-la ANGULO=der_exp_min(1,1)-0.01; for x=1:abs(der_exp_min(1,1)-der_exp_min(aux1,1))/0.01+1 ANGULO=ANGULO+0.01; matriz_der_exp_min(x,1)=ANGULO; matriz_der_exp_min(x,2)=funcao_der_exp_min(ANGULO); if x==1 angulo_fim_comb(2)=matriz_der_exp_min(1,2); else if matriz_der_exp_min(x,2)<angulo_fim_comb(2) angulo_fim_comb(1)=matriz_der_exp_min(x,1); angulo_fim_comb(2)=matriz_der_exp_min(x,2); end; end; end;
angulo_final_da_combustao=angulo_fim_comb(1)
%% Calculando o valor de R² a na expansão isentrópica (após fim da
combustão)
za=round(angulo_fim_comb(1)/intervalo_angulo);
inicio_expansao_isent_1=za*intervalo_angulo; ANGULO=inicio_expansao_isent_1-intervalo_angulo;
% Formando vetor com os dados nos intervalos após fim da expansão
isentrópica determinado em DADOS AUXILIARES for x=1:abs(final_expansao_isent-
inicio_expansao_isent_1)/intervalo_angulo+1
111
ANGULO=ANGULO+intervalo_angulo; LINHA=1+(ANGULO+360)/intervalo_angulo;
lnV_exp(x,1)=ANGULO; lnV_exp(x,2)=lnV(LINHA,2); lnP_exp(x,1)=ANGULO; lnP_exp(x,2)=lnP(LINHA,2); end;
%% Ajuste linear de uma reta na expansão funcao_reta_exp=fit(lnV_exp(:,2),lnP_exp(:,2),'poly1');
% %% Criando matriz com os valores de 0º a 180º % % ANGULO=0-intervalo_angulo; % for x=1:(180/intervalo_angulo+1) % % ANGULO=ANGULO+intervalo_angulo; % LINHA=1+(ANGULO+360)/intervalo_angulo; % % exp_reta_matriz(x,1)=lnV(LINHA,2); % exp_reta_matriz(x,2)=funcao_reta_exp(lnV(LINHA,2)); % exp_reta_matriz(x,3)=ANGULO; % end;
%% calculando o R^2 da aproximação linear
%calculo do valor médio de lnP soma_lnP=0; ANGULO=inicio_expansao_isent_1-intervalo_angulo;
for x=1:abs(final_expansao_isent-
inicio_expansao_isent_1)/intervalo_angulo+1
ANGULO=ANGULO+intervalo_angulo; LINHA=1+(ANGULO+360)/intervalo_angulo;
soma_lnP=soma_lnP+lnP(LINHA,2); end; media_lnP=soma_lnP/x;
% Calculo da soma total do resíduo quadrado e soma total dos quadrados ANGULO=inicio_expansao_isent_1-intervalo_angulo; soma_total_dos_quadrados_exp=0; soma_do_quadrado_dos_residuos_exp=0;
for x=1:abs(final_expansao_isent-
inicio_expansao_isent_1)/intervalo_angulo+1
ANGULO=ANGULO+intervalo_angulo; LINHA=1+(ANGULO+360)/intervalo_angulo;
exp_reta_residuo(x,1)=ANGULO; exp_reta_residuo(x,2)=lnP(LINHA,2)-funcao_reta_exp(lnV(LINHA,2));
soma_do_quadrado_dos_residuos_exp=soma_do_quadrado_dos_residuos_exp+(e
xp_reta_residuo(x,2))^2;
112
soma_total_dos_quadrados_exp=soma_total_dos_quadrados_exp+(lnP(LINHA,2
)-media_lnP)^2;
end;
r2_exp=1-
soma_do_quadrado_dos_residuos_exp/soma_total_dos_quadrados_exp;
% Plotando % plot(exp_reta_matriz(:,1),exp_reta_matriz(:,2)); % hold off
%mostrando r2_comp com 10 dígitos vpa(r2_exp,10)
if size(lnP_exp,1)>pontos_min_para_ajuste_linear_exp if r2_exp>r2_minimo_exp
%% Metodo das DERIVADAS do erro %% compressao
% matriz com a diferença de cada ponto entre o valor experimental de
lnP e o valor de lnP ajustado pela reta
ANGULO=inicio_compressao_isent-intervalo_angulo; for x=1:abs(inicio_compressao_isent-
(final_compressao_isent_1+angulo_a_mais_comp))/intervalo_angulo
ANGULO=ANGULO+intervalo_angulo; LINHA=1+(ANGULO+360)/intervalo_angulo;
erro_comp(x,1)=ANGULO; erro_comp(x,2)=(lnP(LINHA,2)-funcao_reta_comp(lnV(LINHA,2)));
end;
% Plotando % figure; % plot(erro_comp(:,1),erro_comp(:,2))
%% Calculando a derivada do erro_comp. São criadas duas matrizes: 1ª
número de pontos igual ao erro_comp e a 2ª apenas com os pontos
responsáveis por gerar o ajuste linear com critério definido de um R²
mínimo ANGULO=inicio_compressao_isent; for g=1:2;
if g==1 H=final_compressao_isent_1+angulo_a_mais_comp-
intervalo_angulo; else H=final_compressao_isent_1-intervalo_angulo; end;
113
for x=1:abs(inicio_compressao_isent-H)/intervalo_angulo
if g==1 % Calculo aproximado da derivada de erro_comp. Sendo o calculo
feito entre cada dois pontos. O calculo é feito para o mesmo intervalo
de erro_comp. der_comp(x,1)=(erro_comp(x+1,1)-
erro_comp(x,1))/2+erro_comp(x,1); der_comp(x,2)=(erro_comp(x+1,2)-
erro_comp(x,2))/(erro_comp(x+1,1)-erro_comp(x,1)); end;
if g==2 %Calculo aproximado da derivada do erro_comp entre cada dois
pontos para o intervalo de angulos em que se faz o ajuste linear para
a compressão (R²>r2_minimo_comp) der_comp_1(x,1)=(erro_comp(x+1,1)-
erro_comp(x,1))/2+erro_comp(x,1); der_comp_1(x,2)=(erro_comp(x+1,2)-
erro_comp(x,2))/(erro_comp(x+1,1)-erro_comp(x,1)); g2=x; end; end; end;
%% encontrando o maximo valor absoluto da der_comp_1, ou seja, para a
derivada do intervalo de angulo que se conseguiu fazer um ajuste
linear com R²>r2_minimo_comp max_abs_der_comp=abs(der_comp_1(1,2)); for y=2:g2 if abs(der_comp_1(y,2))>max_abs_der_comp max_abs_der_comp=abs(der_comp_1(y,2)); end; end;
%% Plotando a der_comp para um intervalo de ângulo maior que o do
ajuste linear da compressão (R²>r2_minimo_comp) % figure; % plot(der_comp(:,1),der_comp(:,2)) % hold all % %A partir do valor máximo absoluto da função dif_der_comp, plota-se
duas % %retas paralelas ao eixo y, uma com esse valor positivo e outra com
o valor % %negativo % plot(der_comp(:,1),max_abs_der_comp); % plot(der_comp(:,1),-max_abs_der_comp); % hold off
%% encontrando o angulo de inicio da combustão
%diferença de cada ponto entre o valor experimental de lnP e o valor
de lnP ajustado pela reta
% O critério utilizado foi o ângulo em que a der_comp se tornou 10%
maior que o maior der_comp encontrada no intervalo % em que se conseguiu fazer o ajuste linear da compressão
(R²>r2_minimo_comp)
114
y=2; b=0; length_comp=size(der_comp,1); while abs(b)<1.5*max_abs_der_comp y=y+1; if y<=length_comp-2 ang=der_comp(y,1); b=der_comp(y,2); % Como a turbulência gerada pela injeção de combustível afeta o
valor de erro_comp, fazendo este ultimo oscilar bastante, então, para
não encontrar erroneamente o início da combustão será posta mais uma
restrição. A restrição é que a média de alguns pontos de erro_comp
deve ser maior que o valor do maior erro_comp no intervalo usado para
ajustar a reta. Como o salto dado pelo início da combustão é muito
grande e os seguintes também, a restrição a seguir será usada if abs(b)>1.5*max_abs_der_comp if (der_comp(y-
1,2)+der_comp(y,2)+der_comp(y+1,2)+der_comp(y+2,2))/4<2.5*max_abs_der_
comp b=max_abs_der_comp/2; end; end;
else b=999*max_abs_der_comp; end; end;
if y<=length_comp-2 angulo_inicio_da_combustao=ang-intervalo_angulo/2
%% IMEP
%% IMEP _ Trapezoidal rule
for y=1:4 IMEP_matriz(y)=0; for x=(y-1)*180/intervalo_angulo+1:y*180/intervalo_angulo
IMEP_matriz(y)=IMEP_matriz(y)+0.5*(P(x,1)+P(x+1,1))*abs(v(x,2)-
v(x+1,2)) end; end;
IMEP_net=(IMEP_matriz(3)-IMEP_matriz(2))-(IMEP_matriz(4)-
IMEP_matriz(1));
BMEP
FMEP=IMEP_net-BMEP; % Heywod página 47 e 48
%% Pmax e seu respectivo ângulo
ANGULO=-360-intervalo_angulo; Pmax=P(1); for x=2:(720/intervalo_angulo+1)
115
ANGULO=ANGULO+intervalo_angulo; %angulo de -360 a 360 com
intervalo de intervalo_angulo if P(x)>Pmax Pmax=P(x); ANGmax=ANGULO; end; end; Pmax ANGmax
%% montar uma matriz com os valores calculados do inicio e do final da
combustão para cada curva de pressão que se conseguiu calcular esse
valor pelo método utilizado q=q+1; matriz(q,1)=curva; matriz(q,2)=angulo_inicio_da_combustao;
%degree matriz(q,3)=angulo_final_da_combustao;
%degree matriz(q,4)=IMEP_net*5/10^6;
%kW matriz(q,5)=BMEP;
%kW matriz(q,6)=matriz(q,4)-matriz(q,5);
%kW FMEP matriz(q,7)=BMEP/(((10^5*10^-9*(1200/2)*(1/60)*5)/10^3)*vd);
%bar matriz(q,8)=matriz(q,4)/(((10^5*10^-9*(1200/2)*(1/60)*5)/10^3)*vd);
%bar matriz(q,9)=matriz(q,8)-matriz(q,7);
%bar FMEP matriz(q,10)=Pmax; matriz(q,11)=ANGmax;
%% end; end; end; end;
clearvars -except BMEP curva matriz P q P100 P75 P50 P25 matriz_P25
matriz_P50 matriz_P75 matriz_P100 carga size_xlsx dados_validos_P25
dados_validos_P50 dados_validos_P75 dados_validos_P100 P_IMEP
intervalo_angulo media_desvPad_matriz_P25 media_desvPad_matriz_P50
media_desvPad_matriz_P75 media_desvPad_matriz_P100
end;
%% Exportando os dados da combustão para um arquido do Excel if q>0; if carga==1 matriz_P25=matriz; xlswrite('matriz angulos cada curva',matriz_P25,'P25'); %Exportando
os dados encontrados para um arquivo do excel end;
if carga==2 matriz_P50=matriz; xlswrite('matriz angulos cada curva',matriz_P50,'P50'); %Exportando
os dados encontrados para um arquivo do excel
116
end;
if carga==3 matriz_P75=matriz; xlswrite('matriz angulos cada curva',matriz_P75,'P75'); %Exportando
os dados encontrados para um arquivo do excel end;
if carga==4 matriz_P100=matriz; xlswrite('matriz angulos cada curva',matriz_P100,'P100'); %Exportando
os dados encontrados para um arquivo do excel end; end;
clear matriz;
%% FASE DOIS --> ESCOLHENDO OS DADOS "BONS"
%% Escolhendo a matriz a se trabalhar
varlist={'matriz_P','n','media_matriz_P', 'standard_deviation',
'chauvenet_criterion', 'chauvenet_solve', 'chauvenet_inicio_comb',
'chauvenet_final_comb'}; clear(varlist{:});
if carga==1 matriz_P=matriz_P25; end; if carga==2 matriz_P=matriz_P50; end; if carga==3 matriz_P=matriz_P75; end; if carga==4 matriz_P=matriz_P100; end;
%% Calculando a média do incicio da combustão (coluna 1), média do
final da combustão (coluna 2), média do final da combustão pelo
metodo2 (coluna 3)
n=size(matriz_P,1); for y=1:size(matriz_P,2)-1 soma_matriz_P(y)=0; for x=1:n soma_matriz_P(y)=soma_matriz_P(y)+matriz_P(x,y+1); end; media_matriz_P(y)=soma_matriz_P(y)/x;
%% % Calculando a o Standard Deviation standard_deviation(y)=0; for x=1:n
117
standard_deviation(y)=standard_deviation(y)+(matriz_P(x,y+1)-
media_matriz_P(1,y))^2; end; standard_deviation(y)=sqrt(standard_deviation(y)/(n-1));
end;
%% Calculando o valor máximo para (deviation/Satandard_Deviation)
baseando-se no critério de Chauvenet. Critério de Chauvenet pode ser
encontrado no livro: Holman, J. P.; Experimental Methods for
engineers, 6th edition; Página 70) syms chauvenet
chauvenet_solve=solve(0.5*erf(chauvenet/sqrt(2))==(1-
1/(2*n))/2,chauvenet); chauvenet_criterion=vpa(chauvenet_solve);
%% criar uma matriz apenas com os valores que foram considerados
validos simultaneamente no final e no inicio da combustão z=0; for x=1:n chauvenet_inicio_comb=abs((matriz_P(x,2)-
media_matriz_P(1))/standard_deviation(1)); chauvenet_final_comb=abs((matriz_P(x,3)-
media_matriz_P(2))/standard_deviation(2)); chauvenet_IMEP=abs((matriz_P(x,4)-
media_matriz_P(3))/standard_deviation(3));
if chauvenet_inicio_comb<chauvenet_criterion if chauvenet_final_comb<chauvenet_criterion if chauvenet_IMEP<chauvenet_criterion z=z+1; if carga==1 for x1=1:size(matriz_P,2) dados_validos_P25(z,x1)=matriz_P(x,x1); end; end; if carga==2 for x1=1:size(matriz_P,2) dados_validos_P50(z,x1)=matriz_P(x,x1); end; end; if carga==3 for x1=1:size(matriz_P,2) dados_validos_P75(z,x1)=matriz_P(x,x1); end; end; if carga==4 for x1=1:size(matriz_P,2) dados_validos_P100(z,x1)=matriz_P(x,x1); end; end; end; end; end; end;
if z>0; if carga==1
118
n=size(dados_validos_P25,1); for y=1:size(dados_validos_P25,2)-1 soma_matriz_P25(y)=0; for x=1:n
soma_matriz_P25(y)=soma_matriz_P25(y)+dados_validos_P25(x,y+1); end; media_desvPad_matriz_P25(y,1)=soma_matriz_P25(y)/x;
media_desvPad_matriz_P25(y,2)=0; for x=1:n
media_desvPad_matriz_P25(y,2)=media_desvPad_matriz_P25(y,2)+(dados_val
idos_P25(x,y+1)-media_desvPad_matriz_P25(y,1))^2; end;
media_desvPad_matriz_P25(y,2)=sqrt(media_desvPad_matriz_P25(y,2)/(n-
1)); end; media_desvPad_matriz_P25(1,3)=size(matriz_P25,1); media_desvPad_matriz_P25(2,3)=size(dados_validos_P25,1); xlswrite('Resultado.xls',media_desvPad_matriz_P25,'P25'); end;
if carga==2 n=size(dados_validos_P50,1); for y=1:size(dados_validos_P50,2)-1 soma_matriz_P50(y)=0; for x=1:n
soma_matriz_P50(y)=soma_matriz_P50(y)+dados_validos_P50(x,y+1); end; media_desvPad_matriz_P50(y,1)=soma_matriz_P50(y)/x;
media_desvPad_matriz_P50(y,2)=0; for x=1:n
media_desvPad_matriz_P50(y,2)=media_desvPad_matriz_P50(y,2)+(dados_val
idos_P50(x,y+1)-media_desvPad_matriz_P50(y,1))^2; end;
media_desvPad_matriz_P50(y,2)=sqrt(media_desvPad_matriz_P50(y,2)/(n-
1)); end; media_desvPad_matriz_P50(1,3)=size(matriz_P50,1); media_desvPad_matriz_P50(2,3)=size(dados_validos_P50,1); xlswrite('Resultado.xls',media_desvPad_matriz_P50,'P50'); end;
if carga==3 n=size(dados_validos_P75,1); for y=1:size(dados_validos_P75,2)-1 soma_matriz_P75(y)=0; for x=1:n
soma_matriz_P75(y)=soma_matriz_P75(y)+dados_validos_P75(x,y+1); end; media_desvPad_matriz_P75(y,1)=soma_matriz_P75(y)/x;
media_desvPad_matriz_P75(y,2)=0;
119
for x=1:n
media_desvPad_matriz_P75(y,2)=media_desvPad_matriz_P75(y,2)+(dados_val
idos_P75(x,y+1)-media_desvPad_matriz_P75(y,1))^2; end;
media_desvPad_matriz_P75(y,2)=sqrt(media_desvPad_matriz_P75(y,2)/(n-
1)); end; media_desvPad_matriz_P75(1,3)=size(matriz_P75,1); media_desvPad_matriz_P75(2,3)=size(dados_validos_P75,1); xlswrite('Resultado.xls',media_desvPad_matriz_P75,'P75'); end;
if carga==4 n=size(dados_validos_P100,1); for y=1:size(dados_validos_P100,2)-1 soma_matriz_P100(y)=0; for x=1:n
soma_matriz_P100(y)=soma_matriz_P100(y)+dados_validos_P100(x,y+1); end; media_desvPad_matriz_P100(y,1)=soma_matriz_P100(y)/x;
media_desvPad_matriz_P100(y,2)=0; for x=1:n
media_desvPad_matriz_P100(y,2)=media_desvPad_matriz_P100(y,2)+(dados_v
alidos_P100(x,y+1)-media_desvPad_matriz_P100(y,1))^2; end;
media_desvPad_matriz_P100(y,2)=sqrt(media_desvPad_matriz_P100(y,2)/(n-
1)); end; media_desvPad_matriz_P100(1,3)=size(matriz_P100,1); media_desvPad_matriz_P100(2,3)=size(dados_validos_P100,1); xlswrite('Resultado.xls',media_desvPad_matriz_P100,'P100'); end; end;
varlist={'soma_matriz_P25','soma_matriz_P50','soma_matriz_P75','soma_m
atriz_P100'}; clear(varlist{:});
end;
%% Encontrando curva representativa
for carga=1:4
clear dados
if carga==1 P=P25(:,curva); dados=dados_validos_P25; dados_final=media_desvPad_matriz_P25; end; if carga==2 P=P50(:,curva); dados=dados_validos_P50;
120
dados_final=media_desvPad_matriz_P50; end; if carga==3 P=P75(:,curva); dados=dados_validos_P75; dados_final=media_desvPad_matriz_P75; end; if carga==4 P=P100(:,curva); dados=dados_validos_P100; dados_final=media_desvPad_matriz_P100; end;
clear comparacao
for x=1:size(dados,1) comparacao(x,1)=dados(x,1); comparacao(x,2)=(dados(x,10)-dados_final(9,1))/dados(x,10); comparacao(x,3)=(dados(x,11)-dados_final(10,1))/dados(x,11); comparacao(x,4)=(abs(2*comparacao(x,2))+abs(comparacao(x,3)))/3; end;
aux_comparacao(1)=comparacao(1,1); aux_comparacao(2)=comparacao(1,4); for x=2:size(comparacao,1) if comparacao(x,4)<aux_comparacao(2); aux_comparacao(1)=comparacao(x,1); aux_comparacao(2)=comparacao(x,4); end; end;
aux_comparacao
curva_padrao_1(:,1)=P25(:,1); %Angulos
for x=1:4 if carga==x curva_padrao(x)=aux_comparacao(1); curva_padrao_1(:,(x+1))=P(:,1); end; end;
end;
xlswrite('curvas.xls',curva_padrao_1);
121
Anexo III: Calculo da abertura das válvulas de admissão e exaustão
Reta ang 𝑥 𝐷𝑥 𝑇 𝐷𝑒𝑠𝑙(𝐷) [mm]
1 0 25,328 0,000 2,160 0,217
2 5 25,331 0,003 2,166 0,226
3 10 25,339 0,011 2,182 0,247
4 15 25,408 0,080 2,327 0,445
5 20 25,542 0,214 2,607 0,827
6 25 25,823 0,495 3,194 1,632
7 30 26,273 0,944 4,134 2,930
8 35 26,843 1,515 5,326 4,591
9 40 27,437 2,109 6,565 6,338
10 45 28,166 2,838 8,085 8,504
11 50 28,895 3,567 9,608 10,700
12 55 29,560 4,232 10,996 12,724
13 60 30,162 4,834 12,256 14,576
14 65 30,666 5,338 13,313 16,142
15 70 31,083 5,755 14,187 17,445
16 75 31,404 6,076 14,863 18,455
17 80 31,623 6,295 15,324 19,147
18 85 31,721 6,393 15,530 19,458
19 90 31,697 6,369 15,479 19,381
20 95 31,551 6,223 15,172 18,919
21 100 31,286 5,957 14,614 18,082
22 105 30,905 5,576 13,813 16,886
23 110 30,413 5,085 12,781 15,353
24 115 29,765 4,437 11,425 13,352
25 120 29,069 3,740 9,969 11,225
26 125 28,342 3,014 8,452 9,031
27 130 27,549 2,221 6,798 6,668
28 135 26,835 1,507 5,308 4,566
29 140 26,300 0,972 4,192 3,010
30 145 25,929 0,601 3,416 1,938
31 150 25,645 0,317 2,823 1,123
32 155 25,448 0,120 2,411 0,560
33 160 25,333 0,005 2,170 0,231
34 165 25,328 0,000 2,160 0,217
122
Anexo IV: Programa para calcular o deslocamento das válvulas
%% Programa para calcular deslocamento de D
clear clc
deltax=xlsread('deltaX.xls');
a=162.24; b=19.14; c=27.27; l=degtorad(2); m=degtorad(83); syms T
for z=1:size(deltax,1) dx=deltax(z,4);
f=-162.24^2+(162.24*sin(degtorad(83))-
27.27*sin(degtorad(2))+27.27*sin(degtorad(T))-dx)^2+(19.14-27.27*(1-
cos(degtorad(2)))+27.27*(1-cos(degtorad(T))))^2;
df=diff(f,1);
a(1)=5; y=1; x=1; while y>0.001 a(x+1)=a(x)-subs(f,T,a(x))/subs(df,T,a(x)); x=x+1; y=abs(a(x)-a(x-1)); end;
deslD(z,1)=deltax(z,1); deslD(z,2)=deltax(z,2); deslD(z,3)=dx; deslD(z,4)=a(x); deslD(z,5)=43.77*(sin(degtorad(30))-sin(degtorad(30-(a(x)-2)))); deslD(z,6)=43.77*(sin(degtorad(30))-sin(degtorad(30-(a(x)-2))))/0.487;
clearvars a f df
end;
xlswrite('deslD.xls',deslD)
