Pontifícia Universidade Católica de Minas Gerais
Programa de Pós-Graduação em Engenharia Mecânica
________________________________________________
Dissertação de Mestrado
“ESTUDO DE VIABILIDADE DA
RECUPERAÇÃO DE CALOR DOS GASES DE
EXAUSTÃO EM MOTORES PARA
REFRIGERAÇÃO DE CARGAS TÍPICAS EM
MEIOS DE TRANSPORTE RODOVIÁRIO”
VALBERT GARCIA ASSUMPÇÃO
ORIENTADORA: Profª. Elizabeth Marques Duarte Pereira, Drª.
Março de 2004
Pontifícia Universidade Católica de Minas Gerais
Programa de Pós-Graduação em Engenharia Mecânica
________________________________________________
Dissertação de Mestrado
“ESTUDO DE VIABILIDADE DA
RECUPERAÇÃO DE CALOR DOS GASES DE
EXAUSTÃO EM MOTORES PARA
REFRIGERAÇÃO DE CARGAS TÍPICAS EM
MEIOS DE TRANSPORTE RODOVIÁRIO”
VALBERT GARCIA ASSUMPÇÃO
Dissertação apresentada ao Programa de
Pós-graduação em Engenharia Mecânica
como parte dos requisitos para obtenção do
título de Mestre em Ciências em Engenharia
Mecânica.
Banca Examinadora:
Profª. Elizabeth Marques Duarte Pereira – PUC Minas – Presidente, Orientadora
Profª. Andréa Teixeira Charbel, Drª. – UNI - BH – Examinadora Externa
Prof. Luiz Machado, Dr. UFMG – Examinador Externo
Prof. Sérgio de Morais Hanriot, Dr. – PUC Minas – Examinador Interno
Belo Horizonte, 31 de março de 2004
FICHA CATALOGRÁFICA
Elaborada pela Biblioteca da Pontifícia Universidade Católica de Minas Gerais
A851e
Assumpção, Valbert Garcia
Estudo de viabilidade da recuperação de calor dos gases de exaustão em
motores para refrigeração de cargas típicas em meios de transporte rodoviário.
/ Valbert Garcia Assumpção. Belo Horizonte, 2004.
141f. : Il.
Orientadora: Elizabeth Marques Duarte Pereira,
Dissertação (Mestrado) – Pontifícia Universidade Católica de Minas Gerais.
Programa de Pós-Graduação em Engenharia Mecânica.
1. Refrigeração. 2. Absorção. 3. Climatização. 4. Energia – Consumo.
I. Pereira, Elizabeth Marques Duarte. II. Pontifícia Universidade Católica de
Minas Gerais. Programa de Pós-Graduação em Engenharia Mecânica. III. Título.
CDU: 621.56
Dedico este trabalho á minha esposa, Fernanda,
pelos inúmeros momentos compartilhados e vividos;
aos meus pais pelo apoio incontinente e inspirador;
aos meus avós, pelo exemplo e trajetória de vida.
i
AGRADECIMENTOS
Agradeço imensamente à professora Elizabeth Marques Duarte Pereira, minha
orientadora, pelos ensinamentos e pelo convívio, além de todo o apoio recebido durante
o desenvolvimento deste trabalho.
Agradeço a todo o pessoal que conheci no Green-PUC, pelo apoio e pelo aprendizado
compartilhado; à Dulce, Geraldo, Breno, Daniela e Elisiane, dentre tantos outros, e em
especial ao Daniel Teixeira Gervásio e Vinícius Meireles Ciríaco, pela colaboração no
desenvolvimento do programa computacional apresentado neste trabalho. Agradeço aos
professores do Instituto Politécnico, Célia Mara Sales Buonicontro e Carlos Henrique
Guerra Schettino pela oportunidade e colaboração durante meu estágio de docência e ao
pessoal das oficinas do departamento de Engenharia Mecânica, pelo apoio no uso dos
laboratórios. Agradeço ao professor Lauro de Vilhena Brandão Machado Neto pelo
grande apoio e também a professora Julia Maria Garcia Rocha pelas diversas
contribuições e comentários relevantes para o desenvolvimento deste trabalho.
Agradeço ao pessoal da secretaria do programa de pós-graduação, em especial à
Valéria, e ao Jomar, do departamento de informática, por todo o suporte e atenção
dedicados, e também ao pessoal da biblioteca PUC-Minas.
Agradeço a todos aqueles que de maneira direta ou indireta colaboraram para que este
trabalho pudesse ser concretizado. Agradeço ao Ewerton Salles e ainda ao engenheiro
da empresa Recrusul, Sr. José Armindo Reichert pela gentileza e imensa colaboração.
Agradeço também aos professores Ronaldo Darwich Camilo e Eduardo Schirm
(CEFET) e Ramon Molina (UFMG), pelo apoio e pelas contribuições.
Agradeço aos colegas e professores do programa de pós-graduação, que partilharam
juntamente a oportunidade de desenvolvimento que o convívio e o trabalho
proporcionaram, e também aos coordenadores Perrin Smith Neto e José Ricardo Sodré;
e ainda á FAPEMIG, pela bolsa de estudo. Por fim, agradeço a todos os meus amigos e
familiares, pela cooperação, pelo incentivo e pelo carinho.
ii
RESUMO
Este trabalho buscou avaliar a viabilidade de aproveitamento da energia contida nos
gases de exaustão, produzidos em um motor diesel, para obtenção de efeito de
refrigeração dentro de um baú frigorífico. Essa energia pode ser utilizada como aporte
energético em ciclos de refrigeração por absorção, gerando o efeito de refrigeração para
cargas em transporte rodoviário. Desta forma, a potência consumida no funcionamento
do equipamento de refrigeração convencional pode ser economizada, gerando uma
diminuição da quantidade de rejeitos tóxicos lançados na atmosfera e redução dos
custos operacionais durante seu transporte. Um modelo matemático foi desenvolvido
para determinação do efeito de refrigeração requerido em função do tipo de carga,
temperatura de armazenamento, parâmetros construtivos das carrocerias e condições
climáticas típicas. O programa computacional foi implementado, utilizando-se a
plataforma MATLAB. Foram realizados ensaios experimentais em laboratório e em
campo para comparação e validação do modelo proposto. Os resultados mostraram-se
bastante satisfatórios, com desvios absolutos máximos entre os valores das temperaturas
medidas e simuladas da ordem de 2,5 ºC.
Palavras-chave: refrigeração - absorção; modelagem matemática; climatização; uso eficiente de
energia
.
iii
ABSTRACT
This study aims to evaluate the viability of exhaust gases heat recovery, generated by
the combustion in a Diesel motor, in order to produce a cooling effect inside a
refrigeration truck. The heat can be used as energetic input for an absorption
refrigeration cycle to produce the required cooling effect. Thus, the power consumed
during the operation of the conventional refrigeration equipment can be saved, reducing
toxic waste disposal in the atmosphere and the operational costs during the
transportation. A mathematic model has been developed to define the cooling effect
required, related to the load, storage temperature, constructive parameters of the back of
the truck and typical climatic conditions. The software was implemented using
MATLAB framework. Experimental tests were realized at the laboratory and in the
field for comparison and to validate the proposed model. The results were considered
satisfactory, with maximum absolute error between experimental and simulated
temperatures of no more than 2,5ºC.
Key-words: refrigeration - absorption; mathematic modeling; climatization; energy efficient use
iv
SUMÁRIO
CAPÍTULO 1 – INTRODUÇÃO
1.1 – Motivação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.2 – Objetivos gerais e específicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.3 – Estado da Arte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.4 – Escopo da dissertação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
CAPÍTULO 2 – FUNDAMENTAÇÃO SOBRE REFRIGERAÇÃO
2.1 – O ciclo de refrigeração por compressão de vapor . . . . . . . . . . . . . . . . . . . . . .
5
2.2 – O ciclo de refrigeração por absorção . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.3 – O ciclo de absorção contínuo com simples estágio . . . . . . . . . . . . . . . . . . . . . .
9
2.3.1 - O ciclo água-brometo de lítio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2.3.2 – O ciclo amônia-água
.....................................
16
2.4 – Ciclos de absorção a vapor multiestágio e complexo . . . . . . . . . . . . . . . . . . . .
20
2.4.1 – O ciclo com geração por duplo efeito . . . . . . . . . . . . . . . . . . . . . . . .
21
2.4.2 – O ciclo de refrigeração com duplo efeito . . . . . . . . . . . . . . . . . . . . . .
22
2.4.3 – O ciclo com geração em cascata
25
............................
CAPÍTULO 3 – MODELAGEM MATEMÁTICA
3.1 – Recuperação de energia dos gases de exaustão . . . . . . . . . . . . . . . . . . . . . . . .
28
3.1.1 – Caracterização do motor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
v
3.1.2 – Energia nos gases de exaustão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
3.1.3 – Energia disponível ao gerador . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
3.2 – Caracterização do equipamento de absorção . . . . . . . . . . . . . . . . . . . . . . . . . .
34
3.3 – Carga térmica
35
..................................................
3.3.1 – Considerações gerais
.....................................
35
3.3.2 – Desenvolvimento do modelo matemático . . . . . . . . . . . . . . . . . . . . .
37
3.3.3 – Método Numérico
.......................................
49
..................................................
67
3.4 – Decisão final
CAPÍTULO 4 – VALIDAÇÃO DO MODELO E DISCUSSÃO DOS
RESULTADOS EXPERIMENTAIS
4.1 – Materiais e métodos
.............................................
68
4.2 – Convecção natural – ensaio 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
4.3 – Convecção natural – ensaio 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
78
4.4 – Convecção natural + radiação– ensaio 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
80
CAPÍTULO 5 - SIMULAÇÃO MATEMÁTICA - ANÁLISE DE
RESULTADOS
5.1 – Estudo de caso – Resfriamento de Carne . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
85
CAPÍTULO 6 – CONCLUSÕES E CONSIDERAÇÕES FINAIS
98
Referências Bibliográficas
...........................................
100
Anexos
Anexo A.1 – Listagem do programa
....................................
vi
108
Nomenclatura
Variáveis
a, b, c
coeficientes de correlação de Bennett
Aen
área da superfície externa da parede n, m²
Asup
área da superfície, m²
ai ,n
coeficientes para cada componente i dos gases de exaustão
an
coeficientes do polinômio
Bi
número de Biot
cn
concentração do produto na solução no ponto n
cp
calor específico a pressão constante para cada componente, J / kg ⋅ K
c 0p ,i
calor específico em base molar à pressão constante para o estado padrão de
cada componente i , J / kgmol ⋅ K
COP
coeficiente de performance
C pv
capacidade calorífica, kJ / kg ⋅ K
d
dia do ano
d mes
número de dias do mês
E
força eletromotriz, V
E& ent
energia afluente ao volume de controle
vii
En
espessura da parede n, m
∆X
espessura da parede, m
E pp
espaçamento entre a superfícies da carga e as superfícies internas das paredes,
m
E& ac
variação da energia acumulada no volume de controle
g
aceleração da gravidade, m/s²
G SC
constante solar, 1.367 W/m²
h
entalpia específica, kJ / kg ⋅ K
h
coeficiente de transferência de calor por convecção médio, W / m 2 ⋅ K
H
altitude, km
H
radiação global diária em média mensal para superfície horizontal, MJ/m²
H0
radiação extraterrestre incidente sobre uma superfície horizontal, MJ/m²
hr
coeficiente de transferência de calor por radiação, W / m 2 ⋅ K
hs
hora, h
i
número de passos
I
radiação solar em média horária, MJ/m²
k
condutividade térmica, W / m ⋅ K
L
comprimento característico da superfície, m
m&
vazão mássica do fluido, kg/s
n
número de mols, kgmol
n
número de horas de insolação diária em média diária
viii
nmes
número de horas de insolação diária em média mensal
N
número teórico de horas de insolação diária em média mensal
ni
número de moles do componente i , kgmol
Nu
número de Nusselt médio
P
pressão, kPa
p
perímetro da superfície, m
Pr
número de Prandtl
q
fluxo de transferência de calor, W/m²
q&
taxa de transferência de calor específico, kJ/kg
qS
energia proveniente da radiação solar, W
Q&
taxa de transferência de calor, W
R
constante universal dos gases, J / kgmol ⋅ K
Ra
número de Rayleigh
Re
número de Reynolds
r
número de Fourier
rT
distribuição temporal da radiação global em superfícies horizontais
∆t
intervalo de tempo, s
T
temperatura, K
Tc
temperatura da carga, K
Tm−1
temperatura do ar interno, K
ix
Tm ,n
temperatura superficial interna da parede n, K
Tm+1,n
temperatura superficial externa da parede n, K
T∞
temperatura ambiente, °C
Tamb
temperatura ambiente em média horária, K
Tceu
temperatura do céu, K
Tdp
temperatura do ponto de orvalho, °C;
Tf
temperatura média da película do fluído, K
Tgas
temperatura da mistura dos gases de exaustão, K
Ti
temperatura inicial do corpo de prova, K
Tmax .
temperatura máxima diária em média mensal, K
Tmín.
temperatura mínima diária em média mensal, K
Ts
temperatura mínima de saída da mistura dos gases de exaustão, K
Tsup
temperatura da superfície, K
UR
umidade relativa do fluído, decimal
V
velocidade, m/s
v
volume específico, m³/kg
xi
fração molar de cada componente i
W&
potência, W
Z
cota, m
x
Subscritos:
c
carga
ceu
céu
CO2
dióxido de carbono
dp
ponto de orvalho
f
fluído
gas
gases de exaustão
H 2O
água
O2
oxigênio
D
desabsorvedor
e
entrada
e, n
externo, parede n
E
evaporador
EX
exaustão
g
saturação
G
gerador
G1
gerador 1
G2
gerador 2
i, n
interno, parede n
j
posição do nodo
I
isolamento térmico
xi
m + 1, n
superfície externa da parede n
m −1
ar interno
m, n
superfície interna da parede n
P
bomba
RA
refrigeração por absorção
RC
refrigeração por compressão
s
saída
SC
constante solar
sup, e
superfície externa
sup, i
superfície interna
sup
superfície
Sobrescritos
i
instante anterior
i +1
instante atual
Letras Gregas
α
difusidade térmica, m²/s
αT
absortividade da superfície do teto
xii
β
coeficiente de expansão térmica do fluído, K −1
δ
declinação solar, graus
ε
emissividade da superfície
∈
efetividade do trocador de calor
φ
razão de equivalência combustível/ar
ρ
densidade do fluído, kg/m³
σ
constante de Stefan-Boltzmann, 5.6697 ⋅ 10 −8 W / m 2 ⋅ K 4
µ
viscosidade dinâmica, Ns/m²
µ∞
velocidade do fluído escoando paralelamente à superfície, m/s
ν
viscosidade cinemática, m²/s
ω
valor da hora angular, graus
ωs
hora angular correspondente ao pôr-do-sol, graus
γ
diferença relativa aceitável para temperatura, decimal
ξ
latitude, graus
xiii
CAPÍTULO 1
INTRODUÇÃO
1.1- Motivação
Os aspectos tecnológicos relacionados à conservação das cargas, visando melhorar o
projeto de equipamentos, reduzir as perdas do produto por deterioração, diminuir
custos, mantendo a qualidade dos produtos congelados ou resfriados são conhecidos há
muito tempo ( Reinick, 1994).
Os primeiros caminhões e vagões ferroviários com sistema de refrigeração surgiram na
década de 50. Anteriormente, as mercadorias eram transportadas em recipientes com
gelo e protegidas por carrocerias isolantes, que no início eram fabricadas com cortiça,
palha de arroz ou milho, madeira e isopor. Ao final da década de 60, surgiram as
primeiras carrocerias fabricadas em fibra de vidro e poliuretano injetado.
Reinick (1994) avalia que 40% da produção mundial de alimentos, cerca de 1,7 bilhões
de toneladas, necessitam de refrigeração durante o transporte. No Brasil, a norma
"Transporte de Produtos Alimentícios Refrigerados - Procedimentos e Critérios de
Temperatura", está em vigor desde 29 de junho de 2001, tendo como principal requisito
a uniformidade da temperatura requerida.
Esta norma promove a melhoria da qualidade na cadeia do frio (processadores,
estocadores, fabricantes de equipamentos, além de atacadistas e distribuidores),
regulamentando as condições de transporte dos produtos refrigerados e congelados até o
ponto de venda. O atendimento a tais condições é rigorosamente exigido na exportação
de alimentos, representando um aspecto positivo da carteira de exportações brasileira.
1
Cap. 1 – INTRODUÇÃO
2
Diretamente ligados à qualidade do produto, devem ser agregados os custos
operacionais durante seu transporte, incluindo-se a energia gasta para manter a carga
sob as condições especificadas.
Nos equipamentos de refrigeração em uso atualmente pelos meios de transporte, o ciclo
de compressão mecânica de vapor é amplamente utilizado. Neste sistema, a energia
necessária, na forma de trabalho, para acionamento do compressor é retirada do próprio
motor do veículo, reduzindo sua potência útil, ou de um motor adicional acoplado ao
sistema de refrigeração. Desta forma, verifica-se um acréscimo dos custos operacionais,
decorrente do consumo extra de combustível.
Li (1996) apresenta valores típicos do balanço de energia para motores Diesel, abaixo
discriminados:
Trabalho de eixo: 39,20%
Perda térmica pelos gases de exaustão: 33,20%
Calor na água de resfriamento: 13,84%
Calor para lubrificação de óleo: 4,61%
Perdas por radiação : 9,15%
Assim, constata-se que, uma parcela considerável da energia do combustível, cerca de
1/3, é simplesmente rejeitada na forma de calor nos gases de exaustão, inclusive com
impactos indesejáveis ao meio ambiente. Deve-se ressaltar, entretanto, que nem toda
energia perdida pode ser recuperada. Atualmente, a recuperação de apenas 20% da
energia contida nos gases de exaustão tem se mostrado economicamente viável.
Diante do atual e crescente interesse mundial pelas questões ambientais e pela maior
competitividade industrial e mercadológica, faz-se cada vez mais necessário a busca por
novos equacionamentos e soluções que venham a atender a contínua demanda de
energia imposta pela evolução humana.
Dentro desse cenário, o presente trabalho se propõe a desenvolver uma metodologia de
avaliação do aproveitamento do calor contido nos gases de exaustão, produzidos durante
a combustão em um motor Diesel. Esse calor é utilizado como aporte energético em
ciclos de refrigeração por absorção, gerando o efeito de refrigeração para cargas em
transporte rodoviário.
Cap. 1 – INTRODUÇÃO
3
1.2- Objetivos Gerais e Específicos
Este trabalho propõe a formalização de uma metodologia para avaliação da viabilidade
de recuperação do calor contido nos gases de exaustão de motores Diesel, que operam
sob diferentes regimes de carga e velocidade, em ciclos de refrigeração por absorção.
Objetivos específicos :
-
Desenvolver o modelo matemático para determinação do efeito de
refrigeração requerido em transporte rodoviário em função do tipo de carga,
temperatura de armazenamento, parâmetros construtivos das carrocerias e
condições climáticas típicas.
-
Desenvolver um programa computacional para implementação do modelo e
automatização dos cálculos;
-
Validar o modelo desenvolvido para condições controladas em laboratório e
medidas experimentais em campo.
-
Desenvolver uma metodologia de avaliação da viabilidade da recuperação
dos gases de exaustão do motor Diesel em função da carga de refrigeração
requerida no transporte rodoviário de perecíveis.
1.3- Estado da Arte
Horuz (1999) apresenta um trabalho experimental, onde avalia a aplicabilidade dos
sistemas de refrigeração por absorção de vapor nos veículos de transporte rodoviários, a
partir do aproveitamento da energia disponível nos gases de exaustão do motor, como
fonte de energia. Essa alternativa ao uso dos sistemas convencionais de refrigeração por
compressão de vapor praticamente não compromete o desempenho do motor.
Diehl et al (2001) discute as possibilidades e limitações do aproveitamento do calor
remanescente nos gases de exaustão do motor, após sua passagem pelo sistema
catalisador.
Gordon e Mcbride (1994) apresentam um programa computacional para cálculo das
composições de equilíbrio químico, análise matemática e técnicas para obtenção desse
equilíbrio, além do equacionamento para obtenção das propriedades termodinâmicas e
Cap. 1 – INTRODUÇÃO
4
de transporte para a mistura. A partir de tais composições, promove-se o cálculo teórico
das propriedades do sistema e, portanto, a análise termodinâmica dos equipamentos
utilizados.
Radermacher (2001) apresenta uma fundamentação teórica dos sistemas de absorção em
termos dos ciclos ideais de conversão de energia, com descrição das propriedades dos
fluidos de trabalho e uma análise termodinâmica dos principais sistemas de refrigeração
por absorção.
Chinnappa (1992) apresenta uma descrição geral do funcionamento dos sistemas de
refrigeração por absorção de vapor, revisando as inúmeras pesquisas e trabalhos
desenvolvidos sobre o tema.
Srikhirin et al (2001) apresenta uma revisão da literatura sobre a tecnologia de
refrigeração por absorção, com discussão sobre os vários tipos de sistemas existentes,
pesquisas sobre fluidos de trabalho e dos processos de absorção.
1.4- Escopo da dissertação
No Capítulo 2 é apresentada uma revisão dos fundamentos teóricos sobre ciclos de
refrigeração por compressão e por absorção.
O Capítulo 3 trata da modelagem matemática proposta no desenvolvimento do
programa computacional, juntamente com as considerações para sua implementação.
No Capítulo 4 apresenta-se a validação do modelo matemático, com descrição dos
ensaios experimentais, discussão dos resultados e estudos comparativos dos resultados
experimentais com aqueles obtidos nas simulações.
O Capítulo 5 refere-se à aplicação da metodologia de avaliação da recuperação da
energia para um estudo de caso realizado em um frigorífico localizado na cidade de
Campo Belo, Minas Gerais.
As conclusões finais e recomendações para trabalhos futuros são apresentadas no
Capítulo 6. As referências bibliográficas citadas neste texto e consultadas em sua
preparação constam no final deste manuscrito.
CAPÍTULO 2
FUNDAMENTAÇÃO SOBRE REFRIGERAÇÃO
Neste capítulo apresentam-se os conceitos básicos de processos e equipamentos
térmicos de refrigeração por compressão e absorção. Estes conceitos estão disponíveis
na literatura clássica e foram incluídos neste texto para introdução da nomenclatura e
das definições utilizadas em seu desenvolvimento.
2.1- O ciclo de refrigeração por compressão de vapor
O ciclo de refrigeração por compressão de vapor é o sistema mais utilizado para a
geração do efeito de refrigeração em diversas aplicações práticas. A Figura 2.1 mostra
um esquema com seus equipamentos básicos e o diagrama Temperatura – Entropia
Específica para o ciclo teórico correspondente. Neste caso, não são consideradas perdas
de carga em tubulações e equipamentos pertinentes, assim como irreversibilidades
internas, a exceção feita para a válvula de expansão.
O fluido de trabalho é comprimido a partir do consumo de energia elétrica (w34), sendo
o processo 3-4 considerado isoentrópico no ciclo ideal. O vapor superaquecido gerado
(4) dirige-se ao condensador. A transferência de calor (q41) para a vizinhança ocorre a
pressão constante, incluindo uma região de calor sensível, com temperatura decrescente
(4-4´ ), e uma região em que é rejeitado pelo fluido o calor latente de condensação
(4´-1). O estado à saída do condensador é considerado líquido saturado (1), sendo
expandido isoentalpicamente (1-2) até à pressão do evaporador. Neste equipamento,
ocorre uma absorção de calor também à pressão constante (23) pelo fluido de trabalho,
produzindo o efeito de refrigeração no espaço a ser refrigerado.
5
6
Cap. 2 – Fundamentação sobre refrigeração
(a) Componentes do Ciclo
(b) Diagrama T-s
Figura 2.1 - Ciclo de refrigeração por compressão de vapor
A equação geral da Primeira Lei da Termodinâmica para volumes de controle em
regime permanente é:
.
.
.
Q − W = ∑ ms (hs +
s
onde
.
Vs2
V2
+ g ⋅ Z s ) − ∑ me (he + e + g ⋅ Z e )
2
2
e
.
.
(2.1)
.
Q : a taxa de transferência de calor, W : potência de eixo, m : vazão mássica
do fluido, h: entalpia específica, V: velocidade média, g: aceleração da gravidade, Z:
cota. Os subscritos (s) e (e) referem-se às superfícies de controle de entrada e saída dos
fluxos mássicos, respectivamente.
Para tal situação, a equação da continuidade se reduz a:
.
.
∑ me = ∑ ms
e
(2.2)
s
Desprezando-se as variações das energias cinética e potencial e como cada componente
do ciclo possui uma única entrada e saída, a equação 2.1 pode ser rescrita por unidade
de massa na forma:
q& − ω& = hs − he
(2.3)
7
Cap. 2 – Fundamentação sobre refrigeração
onde q& e ω& representam a transferência de calor e trabalho mecânico específicos,
respectivamente.
Assim, as trocas energéticas específicas para os componentes apresentados na Figura
2.1 são calculadas como:
Evaporador: q 23 = (h 3 - h 2 )
(2.4a)
Condensador: q 41 = (h1 - h 4 )
(2.4b)
Compressor: ω 43 = (h 4 - h 3 )
(2.4c)
Válvula de expansão: h 1 = h 2
(2.4d)
Coeficiente de desempenho para ciclos de refrigeração por compressão de vapor
O coeficiente de desempenho (COPRC) é definido pela razão entre o efeito de
refrigeração obtido durante o processo de evaporação do fluido refrigerante e trabalho
fornecido para sua compressão, considerado em módulo. Assim, tem-se:
COPRC =
q 23
ω 34
=
h3 − h2
h4 − h3
(2.5)
Os valores típicos do coeficiente de desempenho de refrigeradores que operam em
ciclos de compressão de vapor são superiores à unidade, com limites comerciais da
ordem de 4,5, segundo Burghardt (1988).
2.2- Ciclo de refrigeração por absorção
O ciclo de refrigeração por absorção é recomendado para aplicações em que há
disponibilidade de energia proveniente de uma fonte térmica, a qual será entregue ao
fluido de trabalho no gerador. A Figura 2.2 mostra um ciclo típico de operação
contínua, cujo fluido refrigerante é uma mistura amônia – água. A área tracejada inclui
os equipamentos que o diferenciam do ciclo de refrigeração por compressão discutido
na seção anterior.
8
Cap. 2 – Fundamentação sobre refrigeração
Figura 2.2 – Princípio do ciclo de absorção de vapor, operação contínua
Fonte: Adaptado de Sonntag (2003)
A mistura amônia–água, denominada solução forte, recebe calor no gerador, produzindo
vapor de amônia a alta pressão. A seguir, dirige-se ao condensador, válvula de expansão
e evaporador, onde gera o efeito de refrigeração desejado. O fluxo remanescente no
gerador, intitulado solução fraca de amônia, flui para o absorvedor para reabsorção do
refrigerante, retornando a mistura, assim, à sua concentração original. Como a solução
fraca de amônia encontra-se a temperatura mais elevada, promove-se uma transferência
de calor para pré-aquecimento da solução forte à entrada do gerador.
Para determinadas aplicações, como é o caso do refrigerador solar, outro modelo básico
de operação pode ser empregado, o ciclo intermitente. Neste caso, as etapas de geração
e absorção ocorrem com defasagem de tempo.
Coeficiente de desempenho para ciclos de refrigeração por absorção
O coeficiente de desempenho (COPRA) é definido pela razão entre o efeito de
refrigeração obtido durante o processo de evaporação do fluido refrigerante ( q E ) e a
energia térmica consumida no gerador ( qG ). Assim, tem-se:
COPRA =
qE
qG
(2.6)
Cap. 2 – Fundamentação sobre refrigeração
9
Propriedades do refrigerante-absorvente
Dentre as propriedades desejáveis para os refrigerantes, destaca-se o alto calor de
vaporização, baixo calor especifico, além de boa estabilidade térmica. Para o
absorvente, destaca-se a estabilidade química, alto ponto de ebulição e baixa capacidade
calorífica, conforme Chinnappa (1992). Para o uso em ciclos contínuos, uma baixa
viscosidade do absorvente garantirá um menor consumo de energia na bomba, sendo
que os pares água-brometo de lítio ( H 2 O − LiBr ) e amônia-água ( NH 3 − H 2 O ) são os
mais comumente utilizados.
Na Tabela. 2.1 apresenta-se uma comparação entre os vários sistemas de refrigeração
por absorção disponíveis, incluindo-se as características operacionais mais relevantes,
de acordo com Srikhirin et al (2001). As células sombreadas correspondem aos ciclos
que permitem atingir temperaturas inferiores a 0 °C dentro das câmaras de refrigeração,
objeto de estudo desse trabalho.
2.3- O ciclo de absorção contínuo com simples estágio
2.3.1- O ciclo água-brometo de lítio
Este equipamento utiliza uma solução de brometo de lítio como absorvente e água como
refrigerante, sendo que seus componentes principais são mostrados na Fig. 2.3.
Fig. 2.3 - Fluxograma para o ciclo água-brometo de lítio simples estágio.
Fonte: Adaptado de Chinnappa(1992)
No trocador de calor, a solução fraca proveniente do gerador (3) transfere calor para a
solução forte proveniente do absorvedor (1).
Sistema
Ciclo com
simples
efeito
Ciclo com
duplo
efeito
(fluxo em
série)
Ciclo com
duplo
efeito
(fluxo em
paralelo)
Ciclo com
triplo
efeito
Nível de Temperatura operacional
pressão
(°C)
Fonte
Refrigeração
de calor
2
80-110
5-10
Fluído de
trabalho
LiBr/água
Capacidade COP
de
refrigeração
(kW)
2,8 - 28
0,5-0,7
Estágio
atual
Observação
Resfriador
com
água
abundante
1.O mais simples e o mais largamente utilizado
2.Utiliza água como refrigerante, com temperatura de
refrigeração ficando acima de 0 °C
3.Pressão do sistema é negativa
4.É necessário utilizar um absorvedor resfriado a água para
evitar a cristalização em altas concentrações
1.É necessária a retificação do refrigerante
2.A solução usada no equipamento é ambientalmente neutra
3.A pressão de operação é tão alta quanto aquela para NH3
4.Não há problemas com cristalização
5.É apropriado para o uso como bomba de calor devido a
grande faixa de atuação
1.É o ciclo de maior performance disponível comercialmente
2.O calor de condensação do primeiro estágio é usado
como calor de entrada para o segundo estágio
2
120-150
<0
Água/NH3
0,85 - 7,0
0,5
Comercial
3
120-150
5-10
LiBr/água
> 280
0,8-1,2
Resfriador
com
água
abundante
<0
Água/NH3
5-10
LiBr/água
2
4
200-230
Unidade
1.O calor liberado no absorvedor do primeiro estágio é usado
experimental
como calor de entrada para o gerador do segundo estágio
N/A
1,4-1,5
Modelo
1.Sistema de controle de alta complexidade
computacional 2.Deverá utilizar fogo direto, pois a temperatura de entrada
e unidade
requerida é bastante alta
experimental 3.Requer maior manutenção, como resultado da elevada
corrosividade causada pela alta temperatura de operação
10
Fonte: Adaptado de Srikhirin et al (2001)
Cap. 2 – Fundamentação sobre refrigeração
Tabela 2.1 – Estudo comparativo entre ciclos de refrigeração por absorção
Sistema
Nível de Temperatura operacional
pressão
(°C)
Fonte
Refrigeração
de calor
3
Baixa
<0
Ciclo com
efeito
parcial
Sistema
2
Com calor
recuperado
Absorção
3
combinada
com ejetor
(Kuhlenschmid t’s)
Chung's and
Chang's
3
Aphornratana's
3
90-180
<0
Fluído de
trabalho
Água/NH3
Água/NH3
Capacidade COP
de
refrigeração
(kW)
N/A
0,2-0,3
N/A
LiBr/água
5-10
2,0
1.Elimina a cristalização no absorvedor pelo aumento da
pressão, através da operação de um ejetor
2.O efeito refrigerante gerado pelo segundo estágio do
gerador é melhor utilizado para funcionar o ejetor que
para produzir efeito de refrigeração
3.O COP esperado é similar ao sistema convencional
Modelo
1.Uma válvula de solução líquida é substituída por um
computacional ejetor acionado pelo líquido
e unidade 2.A taxa de circulação do líquido é reduzida por causa
experimental
do aumento de refrigerante contido na solução, devido a
maior pressão no absorvedor conseguida pelo ejetor
3.Este sistema é apropriado para o uso com refrigerante
com alta densidade devido as características do ejetor
0,9-1,1
Unidade
1.O ejetor localiza-se entre o gerador e o condensador.
experimental
Isto leva o gerador a operar com alta pressão e
temperatura, fazendo com que a temperatura de
entrada seja significantemente aumentada
2.O COP é elevado tão alto quanto ao duplo efeito, com
o aumento do efeito de refrigeração devido ao uso ejetor
3.A taxa de corrosão pode aumentar devido a alta
temperatura de operação
11
Fonte: Adaptado de Srikhirin et al (2001)
LiBr/água
Observação
Modelo
1.Eficiência pobre e complicada
computacional 2.É apropriado em situações onde o calor fornecido
provem de uma fonte de baixa qualidade
0,5-0,7
Modelo
1.O COP deverá ser melhor que o simples efeito em
computacional 10%
Patente
DMETEG/R2 1
DMETEG/R2 2
180-200
Estágio
atual
Cap. 2 – Fundamentação sobre refrigeração
Tabela 2.1 – Estudo comparativo entre ciclos de refrigeração por absorção (cont.)
Sistema
Eames and
Wu's
Nível de Temperatura operacional
pressão
(°C)
Fonte
Refrigeração
de calor
3
200
5
Sistema com
autocirculação
Yazaki
2
80-110
5
Ciclo de
absorção por
difusão
1
140-200
<0
Ciclo com
membrana
osmótica
Ciclo
absorçãocompressão
2
2
LiBr/água
Capacidade COP
de
refrigeração
(kW)
5,0
1,03
Estágio
atual
Observação
Unidade
1.O jato de vapor atua como uma bomba de calor para
experimental recuperar calor do condensado e prover o retorno ao
gerador
2.O ejetor auxilia na redução da pressão no gerador,
portanto a exaustão do ejetor pode ser usada como calor
de entrada
3.O COP é aumentado devido a redução de calor
rejeitado pelo absorvedor
4.Baixa corrosão devido a baixa temperatura de
operação, < 100 °C
LiBr/água
10,0 - 20,0
0,6
Resfriador à 1.Absorvedor resfriado à água é requerido quando o
água
LiBr é utilizado
2.Não é necessário nenhuma bomba mecânica operando
exceto para a água resfriada e para a água refrigerada
Água, NH3/H2 0,050 - 0,3 0,05-0,2 Refrigerador 1.Ciclo de refrigeração operado com calor direto
ou He
doméstico 2.Pode ser operado onde a energia elétrica é escassa
3.Menor manutenção devido a falta de partes móveis
4.A solução de trabalho é ambientalmente utilizável
Patente
1.O sistema é limitado pela tecnologia da membrana
2.Ciclo operado com calor puro, porem o sistema utiliza
bomba, condensador e trocador de calor para a solução
Vários
> 4,5
Modelo
1.A operação do sistema requer um equipamento
computacional mecânico para acionar o compressor
e unidade 2.O circuito de absorção é utilizado na substituição do
experimental condensador e evaporador do tradicional ciclo de
compressão para redução da taxa de compressão, o que
ajuda na redução da potência de entrada na compressão
12
Fonte: Adaptado de Srikhirin et al (2001)
Fluído de
trabalho
Cap. 2 – Fundamentação sobre refrigeração
Tabela 2.1 – Estudo comparativo entre ciclos de refrigeração por absorção (cont.)
13
Cap. 2 – Fundamentação sobre refrigeração
Os balanços de massa no gerador para o fluxo total e para o brometo de lítio são
expressos como:
m& 2 = m& 5 + m& 3
& 2 = c5 m
& 5 + c3m
&3
c2m
(2.7)
(2.8)
onde
m& n - vazão mássica no ponto n ;
cn - concentração do LiBr na solução no ponto n ;
Considerando-se c 5 = 0 e m& 1 = m& 2 = 1 kg/s , tem-se que :
&3=
m
&5 =
m
c2
&2
⋅m
c3
c3 − c 2
&7 =m
&8
=m
c3
(2.9)
(2.10)
A Primeira Lei da Termodinâmica, equações 2.1 e 2.3, aplicada para os diferentes
equipamentos em regime permanente, quando se desprezam as variações de energia
cinética e potencial, conduz a:
Para o trocador de calor:
m& 3 (h3 − h4 ) = m& 2 (h2 − h1 )
h2 =
m& 3
(h3 − h4 ) + h1
m& 2
(2.11)
Para o gerador:
Q& G = m& 3 h3 + m& 5 h5 − m& 2 h2
(2.12)
Q& E = m& 7 ( h8 − h7 ) = m& 5 ( h8 − h6 )
(2.13)
Para o evaporador:
14
Cap. 2 – Fundamentação sobre refrigeração
onde a numeração em subscrito corresponde aos estados apresentados na Figura 2.3 e os
subscritos G e E designam o gerador e evaporador, respectivamente.
As temperaturas máximas nos processos de fornecimento e rejeição de calor são
escolhidas para que não ocorra cristalização quando a solução fraca é resfriada no
trocador de calor.
Para este ciclo básico, Chinnappa (1992) cita como valores de referência para a solução
de H 2 O − LiBr , os valores típicos de pressão, temperatura e concentração apresentados
na tabela 2.2.
Tabela 2.2 - Características operacionais do ciclo H2O-LiBr
Ponto
1
2
3
4
5
Estado
Pressão Temperatura Concentração do líquido
(kPa)
( °C )
(kg LiBr / kg solução)
0,9
30
0,54
0,54
5,6
77
0,60
37
0,60
solução
solução
solução
solução
H2O vapor
superaquecido
H2O líquido
6
5,6
35
saturado
H2O vapor
0,9
5
8
saturado
Fonte: Adaptado de Chinnappa (1992)
Entalpia
(kJ/kg)
-177,0
-102,3
-82,0
-165,0
2633,6
146,5
2510,8
Os valores da entalpia para os pontos 6 e 8 são obtidos em tabelas termodinâmicas. Para
o cálculo da entalpia no ponto 5, Equação 2.12, considera-se que a temperatura de
início de ebulição no gerador é ( t 2' ), onde t 2' ≠ t 2 , e que a temperatura máxima da
solução é t 3 . Para a entalpia da água( h f ) e do vapor saturado( hg ) à temperatura de
condensação( t 6 ), tem-se:
 t '2 + t 3

h 5 = h g + c pv 
− t6 
 2

onde:
c pv = 1,894 kJ / kg ⋅ °C
(2.14)
Cap. 2 – Fundamentação sobre refrigeração
15
Aplicando então, as equações 2.6 e 2.9 a 2.13, obtém-se:
m& 3 =
0,54
= 0,9 kg/s
0,60
m& 5 =
0,6 − 0,54
= 0,1 kg/s
0,6
h2 = 0,9(−82 + 165) − 177 = -102,3 kJ/kg
h 5 = 2565,4 + 1,894(0,5(65 + 77) − 35) = 2633,6 kJ/kg
Q& G = 0,9( −82) + 0,1( 2633,6) − ( −102,3) = 291,9 kW
Q& E = 0,1(2510,8 − 146,5) = 236,4 kW
Assim, obtém-se:
COP =
236,4
= 0,81
291,9
Uma importante observação quanto à operação de resfriadores com bomba é assegurar
que aqueles não operem com taxas de transferência de calor abaixo do mínimo
necessário, garantindo que a temperatura permaneça suficientemente alta para manter
em ebulição a solução fraca dentro do gerador. Neste ciclo básico, esta temperatura é da
ordem de 65 °C.
Nos grandes sistemas, uma bomba mecânica mantém o movimento da solução entre o
gerador e o absorvedor, além de garantir a recirculação dentro dos componentes
individuais.
A Figura 2.4 mostra a variação do coeficiente de performance em função das
temperaturas de fornecimento(gerador) e rejeição(condensador) de calor para os ciclos
água-brometo de lítio e amônia-água, com estágio simples e trocadores de calor com
solução ideal. A temperatura do evaporador é de 4,4 °C.
16
Cap. 2 – Fundamentação sobre refrigeração
Figura 2.4 – Variação do coeficiente de performance com a temperatura
Fonte: Adaptado de Chinnappa (1992)
2.3.2- O ciclo amônia-água
Neste sistema, mostrado na Figura 2.5, a água é o solvente. Como o vapor que entra em
ebulição no gerador não é composto por NH 3 pura, pois contém água, a porcentagem
desta dependerá da respectiva temperatura de geração.
Nas plantas que utilizam este sistema, um condensador de refluxo(retificador) é
conectado ao gerador, para reduzir o vapor de água que entra no condensador a níveis
aceitáveis. Este procedimento é necessário quando a temperatura do gerador excede a
120 °C (Chinnappa, 1992).
O calor é fornecido no gerador e rejeitado para o ar ambiente no absorvedor e no
condensador. O efeito de refrigeração ocorre no evaporador.
Para este sistema e em condições de equilíbrio, podem ser escritas as equações para
massa e energia do vapor que deixam o gerador, o absorvedor e o condensador de
refluxo.
O balanço de massa para gerador-condensador de refluxo pode ser escrito:
Fluxo total:
m& 2 = m& 7 + m& 3
Amônia:
& 2 = c7m
& 7 + c3m
&3
c2m
17
Cap. 2 – Fundamentação sobre refrigeração
Figura 2.5 – Fluxograma para o ciclo amônia-água com estágio simples.
Fonte: Adaptado de Chinnappa (1992)
Considerando que m& 2 = 1 kg/s, e c 7 = 1.0 ,
Obtém-se que:
&3=
m
(1 − c 2 )
(1 − c 3 )
(2.15)
&7 =
m
c2 − c3
1 − c3
(2.16)
O balanço de massa aplicado ao condensador de refluxo, indica que:
m& 5 = m& 6 + m& 7
e para a amônia:
& 5 = c6m
& 6 + c7m
&7
c5m
A solução do sistema de equações mostra que:
&5 =
m
(c 2 − c 3 )(1 − c 6 )
(c 5 − c 6 )(1 − c 3 )
(2.17)
&6=
m
(c 2 − c 3 )(1 − c 5 )
(c 5 − c 6 )(1 − c 3 )
(2.18)
18
Cap. 2 – Fundamentação sobre refrigeração
Aplicando a Primeira Lei, tem-se:
Trocador de calor:
h2 =
m& 3
(h3 − h4 ) + h1
m& 2
(2.19)
Gerador:
Q& G = m& 3 h3 + m& 5 h5 − m& 6 h6 − m& 2 h2
(2.20)
Q& E = m& 7 (h10 − h9 ) = m& 7 (h10 − h8 )
(2.21)
Evaporador:
Para o trabalho na bomba isoentrópica, tem-se que:
(2.22)
& 2 v1 (P2 − P1 )
WP = m
onde:
v = volume específico da solução
A Tabela 2.3 apresenta alguns valores para o volume específico da mistura águaamônia para diferentes concentrações.
Tabela 2.3. - Dados sobre o volume específico da mistura água-amônia
Temperatura
(°C)
15
30
45
60
75
90
105
120
135
0,3
1116
1127
1141
1157
1174
1193
1213
1236
1262
Concentração ( x 10-6 m3/kg)
0,4
0,5
0,6
1159
1207
1261
1174
1225
1285
1191
1248
1311
1209
1269
1340
1230
1295
1373
1253
1325
1415
1279
1358
1463
1308
1397
1516
1342
1444
-
Fonte: Adaptado de Chinnappa (1992)
0,7
1325
1355
1390
1431
1479
1530
1591
1662
-
Cap. 2 – Fundamentação sobre refrigeração
19
Semelhante ao ciclo H 2 O − LiBr , o comportamento deste ciclo é determinado pela
temperatura de evaporação do evaporador ( t10 ), pelas temperaturas nos processos de
rejeição de calor no absorvedor e no condensador e pela temperatura máxima do ciclo.
Aplicando as Equações. 2.15 a 2.23, e considerando como valores típicos de: t1 = 30 °C,
t 3 = t 5 = 94 °C, t 6 = t 7 = 50 °C, t 8 = 35 °C, t10 = -10 °C (Chinnappa, 1992), obtém-se:
m& 3 =
1 − 0,45
= 0,9 kg/s
1 − 0,39
m& 7 =
0,45 − 0,39
= 0,1 kg/s
1 − 0,39
m& 5 = 0,1(
1 − 0,71
) = 0,12 kg/s
0,954 − 0,71
m& 6 = m& 5 − m& 7 = 0,02 kg/s
h2 = 0,9(192 + 64) + (−117) = 113,4 kJ / kg
Q& G = 0,9(192) + 0,12(1496) − 1(113,4) − 0,02(45) = 238 kW
Q& E = 0,1(1433 − 347,5) = 108,6 kW
WP = 0,001197(1350 − 292) = 1,27 kW
COP =
108,6
= 0,45
238 + 1,27
Alguns valores típicos operacionais para este ciclo são mostrados na Tabela 2.4 e foram
obtidos a partir das equações e dados apresentados acima.
Constata-se que o coeficiente de performance para o ciclo amônia-água é menor,
comparado ao ciclo água-brometo de lítio, conforme Chinnappa (1992).
Tal fato é atribuído à menor temperatura do evaporador no ciclo NH 3 - H 2 O , da ordem
de -10 °C, quando comparada ao valor utilizado no ciclo H 2 O - LiBr (5 °C).
Cap. 2 – Fundamentação sobre refrigeração
20
Deve-se considerar ainda a necessidade de uso de um condensador de refluxo, o que
acarreta redução do volume de refrigerante NH 3 que entra no condensador e
consequentemente, no evaporador
Tabela 2.4 - Características operacionais do ciclo NH 3 - H 2 O
Ponto
Estado
Pressao Temperatura Concentração líquido/vapor x/xv Entalpia
(kg NH3 / kg solução)
(kPa)
( °C )
(kJ/kg)
1
líquido
292
30
0,45
-117,0
2
líquido
1350
0,45
3
líquido
1350
94
0,39
192,0
4
líquido
1350
40
0,39
-64,0
5
vapor
1350
94
0,95
1496,0
6
líquido
1350
50
0,71
45,0
NH3 vapor 1350
7
50
1,00
1346,0
8 NH3 líquido 1350
35
347,5
NH
3
vapor
292
-10
1433,0
10
Fonte: Adaptado de Chinnappa (1992)
2.4- Ciclos de absorção a vapor multiestágio e complexo
Os ciclos de absorção onde dois ou mais circuitos estão presentes são chamados ciclos
multiestágios e podem ser classificados como:
Geração por duplo efeito.
Refrigeração com duplo efeito.
Geração em cascata.
Ciclos regenerativos.
Ciclos mistos.
Cap. 2 – Fundamentação sobre refrigeração
21
2.4.1- O ciclo com geração por duplo efeito
O ciclo com geração por duplo efeito é mostrado na Figura 2.6.
Figura 2.6 - Fluxograma para um ciclo água-brometo de lítio com geração
por duplo efeito. Fonte: Adaptado de Chinnappa (1992).
A solução é aquecida a alta temperatura no gerador( G1 ) por uma fonte externa, sendo
este o aporte energético para o ciclo. O vapor produzido em G1 condensa no gerador
G 2 , devido à sua menor temperatura, sendo que a troca de calor faz efervescer mais
refrigerante.
O COP teórico pode ser calculado pelas equações clássicas, assumindo-se que o vapor
deixa o trocador de calor (um gerador, por exemplo) em condições de equilíbrio. Para
estas condições, Whitlow (1966) obteve o COP igual a 1,43. Para o ciclo com simples
estágio, o COP era da ordem de 0,79.
Temperaturas mais elevadas, em torno de 163 °C, são requeridas para operar este ciclo,
comparada aos níveis entre 80 e 85 °C necessários para operar o ciclo simples,
conforme Chinnappa (1992).
Cap. 2 – Fundamentação sobre refrigeração
22
2.4.2- O ciclo de refrigeração com duplo efeito
A operação deste ciclo é descrita a seguir, para o par NH 3 - H 2 O como refrigeranteabsorvente. A Figura 2.7 mostra o circuito básico, sendo que a concentração da solução
no trecho primário (1-2-3-4) é menor que no trecho secundário (11-12-13-14).
Figura 2.7 - Fluxograma para o ciclo amônia-água com duplo efeito
de refrigeração. Fonte: Adaptado de Chinnappa (1992).
O refrigerante é gerado no circuito primário e entra no evaporador passando pelo
condensador de refluxo, condensador e válvula de expansão, sendo que o primeiro
efeito de refrigeração ocorrerá no evaporador.
Quando deixa o evaporador, a NH 3 entra no circuito secundário para ser absorvida pela
solução dentro do reabsorvedor. Esta solução evapora no desabsorvedor, produzindo o
segundo efeito de refrigeração. Em seguida a solução forte retorna ao circuito primário.
Chinnappa (1992) sugere alguns valores para os limites de temperatura de operação, a
saber:
Evaporador: t 9 = 11 °C
Reabsorvedor: t11 e t14 = (30 – 36) °C
23
Cap. 2 – Fundamentação sobre refrigeração
Desabsorvedor: t12 , t13 = (7 – 13) °C
Absorvedor: t 4 , t1 = (40 – 30) °C
Temperatura máxima no gerador: t 3 = 94 °C
Temperatura do vapor à saída do condensador de refluxo: t 6 , t 7 = 50 °C
Temperatura no condensador: t 8 = 35 °C
A tabela 2.5 fornece os valores típicos para propriedades da NH 3 - H 2 O para o ciclo
com duplo efeito de refrigeração.
Aplicando a Eq. 2.11 e considerando-se as faixas de temperatura apresentadas acima,
temos então os seguintes balanços de energia descritos:
Para os trocadores de calor,
h2 =
m& 3
(h3 − h4 ) + h1
m& 2
(2.23)
=113,4 kJ / kg
h12 =
m& 14
( h13 − h14 ) + h11
m& 11
(2.24)
= -117 kJ / kg
De acordo com a Eq. 2.20, para o gerador tem-se:
Q& G = 238 kW
Reescrevendo a Eq. 2.21, para o evaporador:
Q& E = m& 8 (h10 − h8 )
= 110,75 kW
(2.25)
24
Cap. 2 – Fundamentação sobre refrigeração
Para o desabsorvedor,
Q& D = m& 15 h15 + m& 13 h13 − m& 12 h12
(2.26)
= 141,2 kW
Temos então:
COP =
QE + QD
QG
=
(2.27)
110,75 + 141,2
236
=1,06
Tabela 2.5 - Propriedades da NH 3 - H 2 O no ciclo com duplo efeito de refrigeração
Ponto Estado Pressao Temperatura Concentração x/xv Entalpia Vazão
(kg NH3 / kg solução)
t
mássica
ou
(kg NH3 / kg vapor) (kJ / kg) (kg / s)
(kPa)
( °C )
1
líquido 292
30
0,450
-117,0
1
2
líquido 1350
1
3
líquido 1350
94
0,390
192,0
0,90
4
líquido
40
0,390
-64,0
5
vapor
1350
0,954
1496,0
12
6
líquido 1350
50
0,710
45,0
0,02
7
vapor
1350
50
0,998=1
1346,0
0,10
8
líquido 1350
35
1
347,5
0,10
10
vapor
627
11
1
1455,0
0,10
11 líquido 627
30
0,625
-89,0
12 líquido 292
7
0,625
0,78
13 líquido 292
13
0,570
-188,0
0,68
14 líquido 627
36
0,570
-81,0
15
vapor
292
13
1
1310,0
0,10
17
vapor
292
-10
1
1433,0
Fonte: Adaptado de Chinnappa (1992)
O valor do COP encontrado representa uma substancial evolução frente ao COP
encontrado para o ciclo NH 3 - H 2 O simples estágio, enquanto operados dentro das
mesmas condições.
Cap. 2 – Fundamentação sobre refrigeração
25
O circuito (1-2-3-4-5-6-7-8-16-17) é para um ciclo NH 3 - H 2 O simples estágio com os
mesmos limites de temperatura que o ciclo descrito na sessão 2.3.2, sendo seu
coeficiente de performance teórico igual a 0,45.
A efetividade destes ciclos é dependente do par refrigerante-absorvente, sendo que o
água-amônia é adaptável as mais diversas combinações, segundo Chinnappa (1992).
2.4.3-
O ciclo com geração em cascata.
Neste ciclo mostrado na Figura 2.8 existem dois circuitos de solução, onde no circuito
de baixa pressão ( G1 - A1 ) a solução está mais concentrada que no circuito de alta
pressão ( G 2 - A2 ).
O calor é fornecido em ambos os geradores ( G1 e G 2 ), ocorrendo um simples efeito de
refrigeração no evaporador.
Figura 2.8 – Fluxograma para o ciclo amônia-água com geração em cascata.
Fonte: Adaptado de Chinnappa (1992)
O refrigerante produzido no estágio de baixa pressão ( G1 ) é absorvido no absorvedor do
estágio de alta pressão ( A2 ).
26
Cap. 2 – Fundamentação sobre refrigeração
O COP deste ciclo tem sido calculado utilizando o par NH 3 - H 2 O , para as duas
temperaturas do evaporador, ou seja: -10 °C (para um ciclo de refrigeração) e 5 °C
( para um ciclo de condicionamento de ar).
Na Tabela 2.6 são sugeridos os limites das temperaturas de operação para o par NH 3 H 2O .
Tabela 2.6 – Condições operacionais típicas para o ciclo com geração em cascata
Temperatura do evaporador(°C)
Faixas de temperaturas do absorvedor (°C)
Temperatura do condensador (°C)
Temperatura máxima do gerador (°C)
QG1 (kW)
QG2 (kW)
QE (kW)
COP
Ciclo para
Ciclo para
refrigeração condicionamento
de ar
t 15
-10
5
(t 1 - t 4 )
30 - 40
30 - 40
(t 8 - t 11 )
30 - 37
30 - 37
t 13
35
35
t 3,t 9
65
57
212,6
245,9
157
219,2
108,6
164,3
0,29
0,35
Fonte: Adaptado de Chinnappa (1992)
Reescrevendo as Eqs. 2.19 a 2.21, para os balanços de energia nos trocadores de calor,
geradores e evaporadores, tem-se:
Q& G1 = m& 3 h3 + m& 5 h5 − m& 2 h2 − m& 6 h6
Q& G 2 = m& 10 h10 + m& 12 h12 − m& 9 h9
Q& E = m& 13 ( h15 − h13 )
(2.28)
(2.29)
(2.30)
onde os índices utilizados referem-se aos estados mencionados na Figura 2.8. G1 e G2
são os geradores 1 e 2, respectivamente.
Para o coeficiente de performance, tem-se:
COP =
QE
(QG1 + QG 2 )
(2.31)
Valores típicos para este ciclo, são mostrados na Tabela 2.7, para temperatura de –10°C.
27
Cap. 2 – Fundamentação sobre refrigeração
Tabela 2.7 - Ciclo de refrigeração com geração em cascata
(para temperatura de evaporação igual a -10 °C)
Ponto Estado Pressão Temperatura Concentração x/xv
(kg NH3 / kg solução)
(kPa)
(°C)
1
líquido
292
30
0,450
2
0,450
3
627
65
0,390
4
líquido
40
0,390
5
vapor
627
65
0,979
6
líquido
627
30
0,625
7
vapor
627
30
1,000
8
líquido
627
30
0,625
9
54
0,625
10
1350
65
0,565
11 líquido
37
0,565
12
vapor
1350
65
1,000
13 líquido
1350
35
14
292
15
vapor
292
-10
Fonte: Adaptado de Chinnappa (1992)
Entalpia Vazão mássica
(kJ / kg)
(kg / s)
-117,0
1,000
-10,0
1,000
55,0
0,900
-64,0
1438,0
0,106
-89,0
0,006
0,100
-89,0
0,725
22,0
0,725
56,0
0,625
-74,0
0,625
1387,0
0,100
347,5
0,100
347,5
0,100
1433,0
0,100
CAPÍTULO 3
MODELAGEM MATEMÁTICA
Neste capítulo é apresentada a modelagem matemática e o programa computacional
referentes à avaliação proposta. O programa utiliza a plataforma Matlab como ambiente
de trabalho.
A modelagem segue três linhas principais e paralelas (A,B,C) de análise do problema
físico, ilustradas no fluxograma da Figura 3.1, a saber:
A - quantifica a energia a ser recuperada dos gases de exaustão, visando sua
utilização no gerador do ciclo de refrigeração por absorção.
B - caracteriza os requisitos que devem ser atendidos para a manutenção das
condições apropriadas no interior do espaço refrigerado, que determinarão a
capacidade mínima necessária ao equipamento de refrigeração.
C - quantifica os fluxos de energia para o interior da carroceria, tomando como
base a temperatura necessária à conservação dos produtos.
No desenvolvimento do modelo matemático, considera-se que a carga e o baú foram
termicamente pré-condicionados.
3.1- Recuperação de energia dos gases de exaustão
Em um motor de combustão interna, o processo de queima da mistura ar/combustível
produz uma expressiva quantidade de energia, sendo que apenas uma parcela é
convertida em potência, enquanto o restante é descartado para a atmosfera nos gases de
exaustão.
28
29
Cap. 3 – Modelagem matemática
Início
Caracterização
do motor
Temperatura
dos produtos
Cálculos
termodinâmicos
Seleção dos
ciclos de
absorção
úteis
Energia
disponível
ao gerador
1
COP
Efeito de
refrigeração
Carga de
refrigeração
1
N
Efeito de
refrigeração do
equipamento >=
carga térmica de
refrigeração?
S
Definição
final
Fim
Figura 3.1 – Fluxograma geral do programa
30
Cap. 3 – Modelagem matemática
Esta seção apresenta um roteiro de avaliação do potencial energético disponível nos
gases de exaustão de um motor diesel utilizado em veículos de transporte rodoviário.
O fluxograma da Figura 3.2 exemplifica as etapas seguidas nessa avaliação.
Início
Definir dados característicos do motor,
como valores para temperatura e vazão
dos gases de exaustão, e também razão
de equivalência combustível/ar
Calcular o valor da fração molar de
cada componente dos gases de
exaustão do motor
Calcular o valor do calor específico
de cada componente da mistura
dos gases de exaustão do motor
Calcular o valor do calor específico
para a mistura dos gases de
exaustão do motor
Calcular a quantidade de energia
disponível nos gases de exaustão
do motor
Energia disponível
nos
gases
de
exaustão
Fim
Figura 3.2 – Cálculo do calor disponível nos gases de exaustão.
31
Cap. 3 – Modelagem matemática
3.1.1- Caracterização do motor
Para caracterização do motor, é necessário quantificar-se as grandezas que permitirão
uma avaliação do potencial energético proporcionado pelos gases de exaustão de um
dado motor. Tais informações são obtidas acoplando-se este motor a um dinamômetro e
em um analisador de emissões, que permitam o controle e coleta dos dados operacionais
como vazão, consumo e torque, etc., a serem utilizados na análise termodinâmica.
Análise Termodinâmica
O diesel é um derivado do petróleo e constituí-se basicamente por hidrocarbonetos
parafínicos, olefínicos e aromáticos que apresentam cadeia carbônica contendo de 6 a
30 átomos. Neste estudo, o óleo diesel utilizado é equivalente ao hidrocarboneto
duodecano, conforme Van Wylen et al. (1998), sendo designado pela fórmula C12 H 26 .
Em seguida, é definido o número de mols para cada um dos componentes da mistura
constituinte dos gases de exaustão, cujos principais produtos da combustão são CO2 ,
H 2 O , O2 e N 2 (Heywood, 1988). O número de mols de cada um dos componentes dos
gases de exaustão, é obtido previamente através de um balanço atômico, que inclui a
razão de equivalência combustível/ar do motor.
Para os componentes dióxido de carbono e água, os valores obtidos são
respectivamente:
nCO2 = 12
n H 2O = 13
Para o oxigênio, o número de mols é definido pela Equação 3.1.
nO2 =
[18,5 ⋅ (1 − φ )]
φ
onde:
φ - razão de equivalência combustível/ar;
Para o nitrogênio, o número de mols é definido pela Equação 3.2.
(3.1)
32
Cap. 3 – Modelagem matemática
nO2 =
69,8
(3.2)
φ
O número de mols total da mistura é obtido pelo somatório dos valores ni
correspondente a cada componente i , ou seja:
n = ∑ ni
(3.3)
onde:
n - número de mols total da mistura, kgmol ;
ni - número de mols de cada componente i da mistura, expresso em kgmol ;
Portanto, os valores para a fração molar de cada componente são definidos como:
xi =
ni
n
(3.4)
Cálculo do calor específico
Para as faixas de temperatura e pressão consideradas, a mistura dos gases de exaustão é
modelada como gás ideal.
Para a faixa de temperaturas situadas no intervalo entre 200 K e 1000 K, o calor
específico de cada componente dos gases de exaustão é encontrado pela aplicação da
função polinomial definida na Equação 3.5, proposta por Gordon e McBride (1994).
Para cada componente i no estado padrão à temperatura T(K), o calor específico em
base molar c 0p ,i é definido por:
(
c 0p ,i = R ⋅ ai1 ⋅ T −2 + ai 2 ⋅ T −1 + ai 3 + ai 4 ⋅ T + ai 5 ⋅ T 2 + ai 6 ⋅ T 3 + ai 7 ⋅ T 4
onde:
R
- constante universal dos gases, J / kgmol ⋅ K ;
T
- temperatura dos gases de combustão, K;
ai ,n - coeficientes para cada componente i dos gases de exaustão;
)
(3.5)
33
Cap. 3 – Modelagem matemática
Obtido o valor do calor específico para cada componente, calcula-se o valor dessa
propriedade para a mistura, aqui definida pela Equação 3.6.
(
c p = ∑ xi ⋅ c 0p ,i
)
(3.6)
onde:
c 0p ,i - calor específico em base molar à pressão constante para o estado padrão de
cada componente i , J / kgmol ⋅ K ;
A Tabela 3.1 mostra os valores dos coeficientes ai ,n referentes a cada um dos
componentes i , aplicáveis à faixa de temperatura compreendida entre 200 K e 1000 K.
Tabela 3.1 – Coeficientes para cada componente i dos gases de exaustão.
Produto
CO2
H2 O
N2
O2
ai1
4,943651E+04
-3,947961E+04
2,210371E+04
-3,425563E+04
ai2
-6,264116E+02
5,755731E+02
-3,818462E+02
4,847001E+02
ai3
5,301725E+00
9,317827E-01
6,082738E+00
1,119011E+00
ai4
2,503814E-03
7,222713E-03
-8,530914E-03
4,293889E-03
ai5
-2,127309E-07
-7,342557E-06
1,384646E-05
-6,836301E-07
ai6
ai7
-7,689989E-10
4,955043E-09
-9,625794E-09
-2,023373E-09
2,849678E-13
-1,336933E-12
2,519706E-12
1,039040E-12
Fonte: Adaptado de Gordon e McBride (2002)
3.1.2- Energia nos gases de exaustão
A temperatura destes gases após deixarem o cilindro está situada na faixa entre 200°C
até 500 °C, conforme Heywood (1988). Entretanto, para o aproveitamento da energia
remanescente da queima do combustível, deve-se considerar a temperatura mínima na
qual estes gases circulam dentro do sistema de exaustão. Desse modo, minimiza-se a
formação de componentes corrosivos que possam ocorrer antes que os gases sejam
expelidos com segurança para a atmosfera.
Pela Primeira Lei da Termodinâmica, tem-se:
Q& EX = m& ⋅ c p ⋅ (Tgas − Ts )
(3.7)
34
Cap. 3 – Modelagem matemática
onde:
Q& EX - taxa do fluxo de calor contido nos gases de exaustão, W;
m&
- vazão mássica da mistura dos gases de exaustão, kg/s;
cp
- calor específico à pressão constante para a mistura, J/kg ⋅ K ;
Tgas - temperatura da mistura dos gases de exaustão, K;
Ts
- temperatura mínima de saída da mistura dos gases de exaustão, K;
3.1.3- Energia disponível ao gerador
A transferência da energia contida nos gases de exaustão é feita através de um trocador
de calor, normalmente construído pela inclusão de um conjunto de aletas posicionadas
paralelamente entre si e fixadas transversalmente ao longo da superfície externa do
gerador do equipamento de refrigeração por absorção. O conjunto formado pelas aletas
unidas ao gerador deve estar inserido na região onde o fluxo dos gases de exaustão é
forçado passar.
A taxa efetiva do fluxo de calor entregue ao gerador é dada por:
Q& G = Q& EX ⋅ ∈
(3.8)
onde:
∈ - efetividade do trocador de calor;
3.2- Caracterização do equipamento de absorção
Para caracterizar o equipamento de absorção adequado à manutenção das condições
climáticas no interior da carroceria é necessário definir algumas variáveis, como a
temperatura na qual os produtos devem ser mantidos sob refrigeração, pois ela implica
na definição da capacidade do equipamento de refrigeração e na quantificação da
energia líquida que penetra na carroceria. Essa variável determina a configuração
mínima necessária ao equipamento de refrigeração que atenda a demanda de
refrigeração dentro da carroceria e pode ser conferida na Tabela 2.1.
35
Cap. 3 – Modelagem matemática
& ), aplica-se a definição do coeficiente de
Para obtenção do efeito de refrigeração ( Q
E
performance e obtém-se:
Q& E = COP ⋅ Q& G
(3.9)
3.3- Carga térmica
Na primeira etapa do programa são efetuados os cálculos para obtenção dos valores das
temperaturas do ar ambiente externo e do céu nas condições locais. Para tal, utiliza-se
como dados de entrada, as médias mensais de temperaturas máxima e mínima e
umidade relativa, disponíveis nas Normais Climatológicas, (DNMET, 1992).
O cálculo das propriedades físicas dos fluidos e materiais envolvidos, bem como dos
coeficientes de transferência de calor, constituem a segunda etapa do programa.
A terceira etapa consiste na resolução do sistema de equações algébricas geradas na
modelagem matemática, utilizando o método das diferenças finitas. Os resultados finais
obtidos são as temperaturas do ar interno, da superfície da carga e das superfícies
interna e externa das paredes da carroceria.
A partir dessas temperaturas, efetua-se o cálculo da energia líquida que penetra no
espaço refrigerado, que corresponde à carga de refrigeração a ser promovida.
3.3.1- Considerações Gerais
Geometria
A carroceria para transporte da carga é modelada como um paralelepípedo, cujas
dimensões internas e externas são dados de entrada do programa. Como referência de
orientação à numeração das paredes, considera-se que o observador está posicionado
frontalmente à face externa da parede traseira. A numeração, mostrada na Figura 3.3,
será utilizada posteriormente.
A carroceria é constituída por três pares de paredes paralelas, onde o teto e o piso são
designados pelos índices 1 e 2, respectivamente. Situadas verticalmente, a parede
traseira e a parede frontal, são designadas pelos índices 3 e 4, nesta ordem. As paredes
laterais esquerda e direita são respectivamente designadas pelos índices 5 e 6.
Cap. 3 – Modelagem matemática
36
Figura 3.3 – Esquema da geometria da carroceria
Ambiente Externo
O ambiente no qual a carroceria está inserida é a atmosfera terrestre. Portanto, as
condições de contorno das faces externas das paredes são influenciadas pela própria
atmosfera e pelas trocas radiativas com o céu e a radiação solar incidente. Sua
contribuição pode ser bastante significativa, mas por simplicidade sua influência estará
restringida ao teto da carroceria. Não foram avaliados a ocorrência de dias chuvosos.
Paredes
Construtivamente, as paredes são formadas por três lâminas justapostas. A lâmina
externa é metálica, enquanto a lâmina interna é confeccionada em PRFV(plástico
reforçado com fibra de vidro), possuindo ambas função estrutural. Estas são separadas
por uma terceira lâmina que possui função de isolante e é constituída por espuma rígida
de poliuretano com espessura variável.
A parede traseira é composta por duas portas, que permitem o carregamento e o acesso
ao interior da carroceria. Não foram considerados efeitos das possíveis infiltrações de ar
como decorrência de falhas no sistema de vedação das portas.
37
Cap. 3 – Modelagem matemática
Ambiente interno
O ambiente interno da carroceria contém inicialmente apenas ar atmosférico, que é
resfriado até a temperatura de operação estabelecida, conforme citado anteriormente.
O produto a ser refrigerado, quando for inserido no espaço, deve ser disposto
preferencialmente no centro do mesmo, respeitando um espaçamento mínimo entre a
superfície das paredes internas e a superfície da carga, com o objetivo de favorecer a
circulação do ar internamente.
3.3.2- Desenvolvimento do modelo matemático
O modelo matemático determinará a evolução temporal das temperaturas de cada nodo
avaliado na região de interesse. Este comportamento é fortemente influenciado pelas
condições climáticas externas.
Os valores das temperaturas no primeiro instante de tempo são inicializados como dados
de entrada. Os valores obtidos ao final da primeira simulação, correspondem aos dados
iniciais da simulação seguinte e, assim, sucessivamente.
Os valores atualizados são utilizados também para determinação das propriedades e
coeficientes referentes à simulação seguinte.
Caracterização climática
Inicialmente, é necessário a obtenção dos dados relativos às condições climáticas do
ambiente externo à carroceria. Para tanto, é necessário o cálculo das temperaturas do ar
ambiente externo e do céu, correspondentes ao horário, época do ano e localidade na
qual transcorre a simulação. Na Figura 3.4 é mostrado o fluxograma para o cálculo da
temperatura ambiente.
Para o cálculo da temperatura ambiente em média horária ( T∞ ) foi utilizada a equação
proposta em ASHRAE (1994) e citada por Pereira (2002).
 T − Tmin   Tmax − Tmin 
15 ⋅ [(hs ) − 14] ⋅ π 
T∞ = Tmax −  max
+
 ⋅ cos

2
2
180



 

(3.10)
38
Cap. 3 – Modelagem matemática
onde:
Tmáx. - temperatura máxima diária em média mensal, K;
Tmín. - temperatura mínima diária em média mensal, K;
hs - hora do dia;
Início
Definir a localidade, mês
do ano e o horário inicial
Fornecer valores tabelados
para temperaturas máxima,
mínima e hora
Calcular a variação da
temperatura ambiente para
cada intervalo de tempo
Temperatura
ambiente
Figura 3.4 – Fluxograma para o cálculo da temperatura ambiente
Temperatura do céu
O roteiro para o cálculo da temperatura do céu é mostrado no fluxograma da Figura 3.5.
A pressão de saturação é calculada através da Equação 3.11, obtida a partir do software
Catt2 (1996), para a faixa de temperatura compreendida entre 5°C e 40°C.
Pg = 7 ⋅10 −8 ⋅ T∞3 + 7 ⋅10 −8 ⋅ T∞2 + 6 ⋅10 −5 ⋅ T∞ + 5 ⋅10 −4
(3.11)
39
Cap. 3 – Modelagem matemática
Calcular a pressão de saturação
para o vapor d’água
Definir o valor para
umidade relativa
Calcular a pressão de vapor
para a pressão de saturação
Calcular a temperatura
do ponto de orvalho para
a pressão de vapor
Calcular a
temperatura
do céu
Temperatura
do céu
Figura 3.5 – Fluxograma para o cálculo da temperatura do céu
A pressão de saturação ( Pg ) é, então, utilizada para o cálculo da pressão de vapor ( Pv ),
definida como:
Pv = Pg ⋅UR
(3.12)
onde:
UR - umidade relativa do fluído, decimal;
A temperatura do ponto de orvalho ( Tdp ), definida pela Equação 3.13, foi obtida a partir
do software Catt2 para a faixa de temperatura compreendida entre 5°C e 40°C.
40
Cap. 3 – Modelagem matemática
Tdp = 0,0065 ⋅ Pv5 − 0,1595 ⋅ Pv4 + 1,5601 ⋅ Pv3 − 7,9701 ⋅ Pv2 + 25,739 ⋅ Pv − 12,192
(3.13)
A temperatura do céu proposta por Bliss (1961) e citada por Smith et al. (1994) é:
Tceu

 Tdp 

= T∞ ⋅ 0.8 + 

 250 
1
4
(3.14)
Coeficientes de transferência de calor
A disposição da superfície de cada parede da carroceria dentro do ambiente e o
comportamento do fluido deslocando em relação às mesmas, proporciona a obtenção de
uma ampla gama de valores quando a análise envolve os coeficientes de transferência
de calor por convecção.
Para a determinação dos coeficientes de transferência de calor por radiação, é
importante conhecer a diferença entre as temperaturas das superfícies que trocam calor e
também a composição do material que as compõe.
Coeficientes de transferência de calor por convecção
Preliminarmente, devem ser efetuados os cálculos das propriedades gerais do fluido, no
caso o ar, tanto externas quanto internas. Estas propriedades são a densidade, o calor
específico à pressão constante, a condutividade térmica, a difusidade térmica, a
viscosidade dinâmica e a viscosidade cinemática.
A seqüência de cálculo destas propriedades é apresentada no fluxograma da Figura 3.6.
Para o cálculo das propriedades do fluido, define-se a temperatura média da película do
fluído ( T f ) sobre a superfície como:
Tf =
onde:
Tsup - temperatura da superfície, K;
Tsup + T∞
2
(3.15)
41
Cap. 3 – Modelagem matemática
Fornecer o valor
para temperatura
da superfície
Calcular a temperatura
da película de fluido
Calcular a densidade do
fluido para temperatura
de película
Calcular o calor específico
do fluido para temperatura
de película
Calcular a condutividade
térmica do fluido para
temperatura de película
.
Calcular a difusidade
térmica do fluido para
temperatura de película
Calcular o valor da
viscosidade
dinâmica
do
fluido
para
temperatura de película
Calcular o valor da
viscosidade cinemática
do
fluido
para
temperatura de película
Figura 3.6 – Fluxograma para o cálculo das propriedades gerais
Para a temperatura da película, calcula-se a massa específica do fluído ( ρ ) à partir da
correlação obtida através do software Catt2, na forma:
42
Cap. 3 – Modelagem matemática
ρ = 351,92 ⋅ (T f )−1, 0017
(3.16)
O cálculo do calor específico (cp) é obtido a partir da correlação empírica apresentada
por Rohsenow e Hartnett (1998), como:
) ] [(
[(
)
c p = 1,03409 + − 0,28488708 ∗ 10 −3 ⋅ T f + 0,7816818 ∗ 10 −6 ⋅ T f
2
]+
(3.17)
[(
)
3
] [(
)
+ − 0,4970786 ∗ 10 −9 ⋅ T f + 0,1077024 ∗ 10 −12 ⋅ T f
4
]
Rohsenow e Hartnett (1998) correlacionam o valor da condutividade térmica ( k ) para
película de fluído na faixa temperatura entre 250 K e 1050 K na forma:
[(
) ]
k = −2,276501 ∗ 10 −3 + 1,2598485 ∗10 −4 ⋅ T f +
[(
) ]+ [(1,73550646 ∗10 )⋅ T ]+
+ [(− 1,066657 ∗10 )⋅ T ]+ [(2,47663035 ∗10 )⋅ T ]
+ − 1,4815235 ∗ 10 −7 ⋅ T f
2
(3.18)
f
4
−13
3
−10
5
−17
f
f
Para a difusidade térmica ( α ), tem-se:
α=
k
(ρ ⋅ c p )
(3.19)
Rohsenow e Hartnett (1998) recomendam para cálculo da viscosidade dinâmica ( µ ) na
faixa de temperatura da película de fluido entre 250°C e 600°C:
[
]
µ = −0,98601 + [(9,080125 ∗ 10 −2 )⋅ T f ] + (− 1,17635575 ∗ 10 −4 )⋅ T f 2 +
(3.20)
[(
+ 1,2349703 ∗ 10
−7
)⋅ T ]+ [(− 5,7971299 ∗10 )⋅ T ]
3
4
−11
f
f
Na seqüência é calculado o valor da viscosidade cinemática (ν), definida pela Equação
3.21.
ν=
µ
ρ
(3.21)
Para definir o comportamento do escoamento, calcula-se o número de Reynolds ( Re ),
parâmetro adimensional que representa a razão entre a força de inércia e a força viscosa
na camada limite fluidodinâmica e está definido como:
43
Cap. 3 – Modelagem matemática
Re =
µ ∞ .L
ν
(3.22)
onde:
L - comprimento característico da superfície, m;
Definido o regime de deslocamento do fluido ao longo de uma superfície plana, duas
condições distintas passam a orientar a seqüência de cálculo, uma aplicada ao
escoamento com velocidade forçada e outra com ventilação natural.
A modelagem desses campos de velocidade utiliza o conceito de camada limite, que
envolve a distribuição de temperaturas ao longo da superfície. Ele delimita uma região
compreendida por uma determinada espessura de fluido por sobre e ao longo de uma
superfície plana, onde os gradientes de temperatura e tensão de cisalhamento são
elevados.
O fluxograma com a seqüência dos cálculos dos parâmetros adimensionais da Figura
3.7 é utilizado para regime com escoamento forçado.
O número de Prandtl (Pr) é definido como a razão entre a difusidade de momento e a
difusidade térmica, na forma:
Pr =
ν
α
(3.23)
O parâmetro que fornece o gradiente de temperatura adimensional na superfície e
determina a intensidade da transferência de calor por convecção é obtido
empiricamente. Para superfícies planas em presença de escoamento forçado e regime
turbulento, conforme Incropera e Dewitt (1998), recomenda-se:
Nu = 0,037 ⋅ Re 4 / 5 ⋅ Pr 1 / 3
onde:
Nu : número de Nusselt médio;
(3.24)
44
Cap. 3 – Modelagem matemática
Fornecer o valor de velocidade
e viscosidade cinemática do
fluido para superfície
Definir o valor do comprimento
característico para a superfície
Calcular o valor do número
de Reynolds para a superfície
Calcular o valor do número
de Prandtl para a superfície
Calcular o valor do número
de Nusselt para a superfície
Figura 3.7 – Fluxograma para o cálculo dos parâmetros adimensionais
Para convecção livre, o cálculo do número de Nusselt utiliza equações distintas quando
a superfície está posicionada no plano horizontal ou no plano vertical. Para sua
obtenção, calcula-se o parâmetro adimensional que representa a razão entre o empuxo e
a força viscosa que atua no fluído, conhecido como número de Rayleigh (Ra):
Ra =
g ⋅ β ⋅ (Tsup − T∞ )⋅ L.3
ν ⋅α
(3.25)
onde:
g - aceleração da gravidade, m/s²;
β - coeficiente de expansão térmica do fluido, K −1 ;
O fluxograma da Figura 3.8 mostra a seqüência do cálculo dos parâmetros
adimensionais para convecção livre.
45
Cap. 3 – Modelagem matemática
Para superfície vertical com temperatura uniforme, ausência de escoamento forçado e ao
longo de todo o intervalo de Ra , o número de Nusselt está definido pela Equação 3.26,
na correlação apresentada por Churchill e Chu (1975), e citada por Incropera e Dewitt
(1998).

0,387 ⋅ Ra 1 / 6
Nu = 0,825 +
9 / 16

1 + (0,492 Pr )
[

8 / 27 

2
]
(3.26)
O número de Nusselt para superfícies horizontais está definido nas correlações
propostas por McAdams (1954) e apresentadas por Incropera e Dewitt (1998).
Fornecer o diferencial entre
temperatura da superfície e
temperatura ambiente.
Definir o valor do comprimento
característico para a superfície.
Fornecer o valor de viscosidade
cinemática, difusidade térmica e
coeficiente de expansão para a
temperatura da superfície.
Calcular o valor do
número de Rayleigh
para a superfície.
Calcular o valor do
número de Nusselt
para a superfície.
Figura 3.8 – Fluxograma para o cálculo dos parâmetros adimensionais
Para superfície voltada para cima com temperatura superior à do fluido ou superfície
voltada para baixo com temperatura inferior à do fluido, a correlação é aplicada para
número de Rayleigh compreendido na faixa de 10 4 ≤ Ra ≤ 10 7 , na forma:
46
Cap. 3 – Modelagem matemática
Nu = 0,54 ⋅ Ra 1 / 4
(3.27)
e para Rayleigh na faixa de 10 7 ≤ Ra ≤ 1011 , como:
Nu = 0,15 ⋅ Ra 1 / 3
(3.28)
Para situações onde 10 5 ≤ Ra ≤ 1010 , superfície voltada para cima com temperatura
superior à do fluido ou superfície voltada para baixo com temperatura inferior à do
fluido, o número de Nusselt está definido por:
Nu = 0,27 ⋅ Ra 1 / 4
(3.29)
O comprimento característico (L), conforme Goldstein et al. (1973), Lloyd e Moran
(1974) e citado em Incropera e Dewitt (1998) é definido como:
L≡
Asup
(3.30)
p
onde:
Asup - área da superfície, m²;
p
- perímetro da superfície, m;
O valor médio do coeficiente de transferência de calor por convecção para a superfície
( h ) e condição em estudo é:
h = Nu ⋅
k
L
(3.31)
Coeficientes de transferência de calor por radiação
O cálculo dos coeficientes de transferência de calor por radiação para a superfície de
cada parede permite quantificar a intensidade com que a energia normalmente é perdida
para o céu. O coeficiente de transferência de calor por radiação está representado pelo
símbolo hr , conforme Incropera e Dewitt (1998) e é definido como:
[
]
hr = ε ⋅ σ ⋅ (Tceu ) + (Tsup ) ⋅ (Tceu + Tsup )
2
2
(3.32)
47
Cap. 3 – Modelagem matemática
onde:
ε
- emissividade da superfície;
σ
- constante de Stefan-Boltzmann, 5.6697 ⋅10 −8 W / m 2 ⋅ K 4 ;
Radiação Solar Incidente
Conforme citado anteriormente, uma parcela da energia proveniente da radiação solar
incidente sobre a superfície da carroceria é contabilizada pelo modelo.
Inicialmente, calcula-se o valor da declinação solar ( δ ) aplicando a correlação proposta
por Cooper (1969) e citada por Duffie e Beckman (1991), conforme mostra a Equação
3.33.

 284 + d 

 365 
δ = 23,45° ⋅ sen 2 ⋅ π ⋅ 

(3.33)
onde:
d - dia do ano;
O valor da hora angular correspondente ao pôr-do-sol ( ω S ) é calculado por Duffie e
Beckman (1991) como:
ω s = ar cos(− tan ξ ⋅ tan δ )
(3.34)
onde:
ξ- latitude, graus;
O número de horas de insolação diária em média diária ( n ) é calculado por:
n=
n mes
d mes
(3.35)
onde:
nmes - número de horas de insolação diária em média mensal, fornecido pelo
INMET(1992)
d mes - número de dias do mês;
48
Cap. 3 – Modelagem matemática
O número teórico de horas de insolação diária em média mensal ( N ) é definido como:
N=
2
⋅ω s
15
(3.36)
Para determinar a distribuição temporal da radiação global em superfícies horizontais,
rt , é preciso encontrar o valor da hora angular, apresentada em Pereira (2002).
ω = 15° ⋅ [(h + 0,5) − 12]
(3.37)
Na seqüência, o valor de rT , proposto por Collares-Pereira e Rabl (1979) e citada em
Duffie e Beckman (1991), é:


(cos ω − cos ω s )
π

rT =
⋅ (a + b ⋅ cos ω ) ⋅ 
24
 sen ω −  π ⋅ ω s ⋅ cos ω s
s
 
180





 

 
(3.38)
Os coeficientes a e b são definidos por:
a = 0,409 − 0,5016 ⋅ sen (ω s − 60 )
(3.39a)
b = 0,6609 − 0,4767 ⋅ sen(ω s − 60 )
(3.39b)
O valor da radiação extraterrestre incidente sobre uma superfície horizontal ( H 0 ) ,
proposto por Duffie e Beckman (1991), é:
H0 =
86400 ⋅ G SC 
 2 ⋅ π ⋅ d 
⋅ 1 + 0,033 ⋅ cos
 ⋅ [(cos ξ ⋅ cos δ ⋅ sen ω s ) + (ω s ⋅ sen ξ ⋅ sen δ )]
π
 365 

(3.40)
onde:
G SC - constante solar, 1.367 W/m²;
Baseado no modelo proposto por Bennett (1965) e citado em Pereira (2002), o cálculo
do valor da radiação global diária em média mensal para superfície horizontal ( H ),
aplicável ao Hemisfério Norte, é:
49
Cap. 3 – Modelagem matemática
n


H = H0 ⋅  a + b ⋅ + c ⋅ H 
N


(3.41)
onde:
H
- altitude, km;
O modelo foi adaptado por Nunes et al. (1976), e validado por Nunes et al. (1976) e
Pereira (1991) para utilização no Hemisfério Sul.
Os coeficientes a , b e c , utilizados na Equação 3.41, são chamados coeficientes
empíricos de correlação de Bennett. Estes coeficientes são tabelados e adaptados para
uso no Hemisfério Sul.
Para obter-se a radiação solar em média horária ( I ), utiliza-se:
I = H ⋅ rT
(3.42)
A energia proveniente da radiação solar incidente e absorvida pela superfície externa da
parede número 1 por unidade de tempo (qs) é definida como:
q S = I ⋅ Ae1 ⋅ α abs
(3.43)
onde:
Ae1 - área da superfície externa da parede 1, m²;
α abs - absortividade da superfície externa do teto
3.3.3- Método Numérico
A modelagem está formatada para trabalhar com condições unidimensionais e em
regime transiente. Portanto, supõe–se que gradientes de temperatura ocorram somente
ao longo de uma única direção coordenada e a transferência de calor ocorra
exclusivamente nesta direção.
Para encontrar a carga térmica determina-se a distribuição de temperaturas para as
superfícies das paredes, piso e teto da carroceria e, ainda, das temperaturas do ar interno
e superfície da carga.
Cap. 3 – Modelagem matemática
50
Método das diferenças finitas
Neste método, a equação diferencial parcial da condução do calor é aproximada por um
conjunto de equações algébricas na temperatura em um certo número de pontos nodais
distribuídos numa região. O primeiro passo é representar a equação diferencial da
condução de calor e suas condições de contorno em um conjunto de equações
algébricas.
Obtém-se a forma em diferenças finitas da equação da condução do calor substituindo
as derivadas parciais da temperatura na equação da condução de calor, pelas formas
equivalentes em diferenças finitas, ou escrevendo-se um balanço de energia num
elemento diferencial de volume.
No programa desenvolvido, utilizou-se do método do balanço de energia aplicado à
conservação de energia no volume de controle em torno da região nodal.
Na análise numérica, fez-se a escolha dos pontos discretos posicionando um nodo para
cada uma das superfícies interna e externa das paredes, um para a superfície da carga e
outro para o ar interno. Considera-se que a parede é isotérmica e, portanto, que o nodo
representa a medida da temperatura média da região.
O posicionamento destes pontos é mostrado no circuito térmico equivalente da Figura
3.9, onde a posição dos nodos ao longo do circuito é a superfície interna, designada pelo
índice subscrito m acompanhado do número relativo à parede.
Figura 3.9 – Circuito térmico equivalente
51
Cap. 3 – Modelagem matemática
À direita, estão os nodos superficiais externos indicados pelo índice subscrito m + 1
seguido pelo número da parede. À esquerda, está situado o nodo ar interno, representado
pelo índice m − 1 e, também, o nodo referente à superfície da carga, com o índice c .
São indicados, também ,os espaçamentos longitudinais entre os nodos, onde E n é a
espessura da parede correspondente e E pp corresponde à metade do valor da distância
existente entre as paredes internas e a superfície da carga.
Além da discretização espacial através da escolha dos pontos nodais, o problema é
dividido discretamente no tempo ( t ). Essa condição é atendida pela adoção do inteiro i ,
sendo:
t = i ⋅ ∆t
(3.44)
onde:
i - número de passos para resolução do sistema;
∆t - intervalo de tempo;
O índice i indica a dependência temporal da temperatura, onde a derivada no tempo é
expressa em termos da diferença entre as temperaturas associadas a um instante novo
( i + 1 ) e ao instante anterior ( i ). Assim, os cálculos são realizados em instantes
sucessivos separados pelo intervalo de tempo ∆t . A temperatura é encontrada para
pontos discretos tanto no espaço quanto para o tempo. A Figura 3.10 mostra essa
discretização para os domínios de x e t .
Figura 3.10 – Domínios de x e t em diferenças finitas
52
Cap. 3 – Modelagem matemática
A solução em diferenças finitas é dependente do instante particular em que as
temperaturas são estimadas pelas aproximações em diferenças finitas das derivadas
temporais. Neste programa, é adotado o método explícito de resolução, onde as
temperaturas são estimadas no instante de tempo anterior ( i ).
Critério de estabilidade
Quando aplicado a condições de regime transiente, o método explícito não apresenta
comportamento incondicionalmente estável. Para atender a exigência de estabilidade,
adotou-se o uso do parâmetro que estabelece os valores máximos permitidos para o
valor do intervalo de tempo ( ∆t ).
O valor adequado deste parâmetro é definido por uma relação limite que envolve o
número de Fourier em diferenças finitas e a reunião de todos os termos relacionados ao
nodo em estudo para obtenção de um coeficiente.
Neste programa definiu-se que o critério é atendido quando o coeficiente associado ao
mesmo nodo é maior ou igual a zero para o instante de tempo anterior.
A forma do número de Fourier em diferenças finitas é definida pela Equação 3.45.
r=
α ⋅ ∆t
(∆X )2
(3.45)
Equações algébricas
A determinação numérica da distribuição de temperatura impõe uma equação de
conservação apropriada para cada um dos pontos nodais de temperatura desconhecida.
O conjunto de equações resultantes é então resolvido sucessivamente na temperatura de
cada nodo.
Os balanços de energia foram formulados supondo-se que todos os fluxos de calor estão
direcionados para dentro do nodo. Fisicamente, esta concepção não ocorre. Entretanto,
se as equações das taxas são expressas coerentemente com esta hipótese, são obtidas as
formas corretas da equação em diferenças finitas.
53
Cap. 3 – Modelagem matemática
O método do balanço de energia expresso na forma geral, quando aplicado a um volume
de controle em torno de um nodo e considerando-se a variação da energia térmica
acumulada, é:
E& ent = E& ac
(3.46)
onde:
E& ent - energia efluente ao volume de controle;
E& ac - variação da energia acumulada no volume de controle;
Um corte da seção transversal de uma das paredes da carroceria é mostrado na Figura
3.11, onde as superfícies interna e externa desta parede estão indicadas respectivamente
como A e e A i . Sobre estas superfícies estão situados dois nodos, representados
respectivamente pelos pontos Te e Ti . A espessura da parede está indicada por ∆X .
No desenvolvimento das equações algébricas para os nodos situados em fronteiras
sujeitas a transferência de calor por convecção ou com fluxos de calor determinados, o
balanço de energia é aplicado ao elemento de volume relativo à metade do valor de
espessura da parede.
Figura 3.11 – Esquema mostrando um corte transversal da parede
Cap. 3 – Modelagem matemática
54
Nodos superficiais externos
Os vários planos representados pelas superfícies externas das paredes da carroceria são
mostrados na Figura 3.12. Situados sobre cada um dos planos estão localizados os
nodos correspondentes a cada uma das seis superfícies.
Figura 3.12 – Planos representativos das superfícies externas das paredes
Para exemplificar, a energia que aflui ao nodo superficial externo da parede 1 provém
do interior da parede por transferência de calor por condução, do ar ambiente externo
através da transferência de calor por convecção e através de transferência de calor por
radiação que ocorre com o céu. O balanço de energia para este nodo considera também
o calor absorvido durante o período de incidência da radiação solar sobre a superfície.
Aplica-se o balanço de energia ao nodo da superfície externa da parede número 1 para
obter-se a forma explícita da equação em diferenças finitas (Özişik, 1990), definida por:
 
 E1 ⋅ hre,1 E1 ⋅ he ,1    i
Tmi ++11,1 = (2 ⋅ r1 ) ⋅ Tmi ,1 + 1 − 2 ⋅ r1 ⋅ 
+
+ 1   ⋅ Tm+1,1 +
kI
 
 kI
  
 2⋅r ⋅ E 
 2 ⋅ r1 ⋅ E1 ⋅ hre,1 
 2⋅r ⋅ E ⋅h 
 ⋅ TCÉU +  1 1 e,1  ⋅ T∞ +  1 1  ⋅ α e ⋅ q
+ 
 k ⋅A 
kI
kI




 I e,1 
(3.47)
Cap. 3 – Modelagem matemática
55
onde:
r1
- critério de estabilidade para o nodo no instante i ;
E1
- espessura da camada de isolamento térmico, m;
hre,1 - coeficiente de transferência de calor por radiação, W / m 2 ⋅ K ;
kI
- condutividade térmica do isolamento térmico, W / m ⋅ K ;
Ae ,1 - área superficial externa, m²;
he,1
- coeficiente de transferência de calor por convecção, W / m 2 ⋅ K ;
Tmi ,1 - temperatura superficial interna, parede 1, K ;
Tmi +1,1 - temperatura superficial externa, parede 1, K;
α e - absortividade superficial externa, parede 1;
qS
- fluxo de calor gerado pela radiação solar incidente, W/m²;
obs.: para índice subscrito, e indica superfície externa e 1 indica parede;
para índice sobrescrito, i indica o instante.
Para o nodo situado na superfície externa da parede número 2, a aplicação do balanço de
energia permite a obtenção da equação em diferenças finitas em sua forma explícita,
conforme mostra a Equação 3.48. Neste nodo a energia contabilizada aflui do interior da
parede por transferência de calor por condução e do ar ambiente externo através da
transferência de calor por convecção.
 
 E ⋅h
 
Tmi ++11, 2 = (2 ⋅ r2 ) ⋅ Tmi , 2 + 1 −  2 ⋅ r2 ⋅  2 e, 2 + 1  ⋅ Tmi +1, 2 +
 
 kI
 
 2 ⋅ r2 ⋅ E 2 ⋅ he, 2 
 ⋅ T∞
+ 
kI


(3.48)
56
Cap. 3 – Modelagem matemática
onde:
r2
- critério de estabilidade para o nodo no instante i ;
E2
- espessura da camada de isolamento térmico, m;
he, 2
- coeficiente de transferência de calor por convecção, W / m 2 ⋅ K ;
Tmi , 2 - temperatura superficial interna, parede 2, K ;
Tmi +1, 2 - temperatura superficial externa, parede 2, K;
Para os nodos localizados sobre as superfícies externas das paredes de número 3, 4, 5 e
6, a forma explícita da equação em diferenças finitas é similar, respeitados os índices
apropriados.
O balanço de energia considera a transferência de calor por convecção com o ar
ambiente externo e a transferência de calor por condução, que ocorre com o interior da
parede. Além destes, contabiliza também a transferência de calor por radiação que
ocorre com o céu. Para o nodo localizado sobre a superfície da parede 3, tem-se:
 
E ⋅h
 E ⋅ hr
  
Tmi ++11,3 = (2 ⋅ r3 ) ⋅ Tmi ,3 + 1 − 2 ⋅ r3 ⋅  3 e,3 + 3 e,3 + 1   ⋅ Tmi +1,3 +
kI
 
 kI
  
 2 ⋅ r3 ⋅ E3 ⋅ hre,3 
 2 ⋅ r3 ⋅ E3 ⋅ he,3 
 ⋅ TCÉU + 
 ⋅ T∞
+ 
kI
kI




(3.49)
Para o nodo da superfície externa da parede 4, tem-se:
 
E ⋅h
 E ⋅ hr
 
Tmi ++11, 4 = (2 ⋅ r4 ) ⋅ Tmi , 4 + 1 − 2 ⋅ r4 ⋅  4 e, 4 + 4 e, 4 + 1  ⋅ Tmi +1, 4 +
kI
 
 kI
 
 2 ⋅ r4 ⋅ E 4 ⋅ hre, 4 
 2 ⋅ r4 ⋅ E 4 ⋅ he, 4
 ⋅ TCÉU + 
+ 
kI
kI




 ⋅ T∞

(3.50)
Cap. 3 – Modelagem matemática
57
O nodo da superfície externa da parede 5 está definido como:
 
E ⋅h
 E ⋅ hr
 
Tmi ++11,5 = (2 ⋅ r5 ) ⋅ Tmi ,5 + 1 − 2 ⋅ r5 ⋅  5 e,5 + 5 e,5 + 1  ⋅ Tmi +1,5 +
kI
 
 kI
 
 2 ⋅ r5 ⋅ E5 ⋅ hre,5 
 2 ⋅ r5 ⋅ E5 ⋅ he,5 
 ⋅ TCÉU + 
 ⋅ T∞
+ 
kI
kI




(3.51)
O nodo da superfície externa da parede 6 é:
 
 E6 ⋅ hre,6 E6 ⋅ he ,6
  
Tmi ++11, 6 = (2 ⋅ r6 ) ⋅ Tmi , 6 + 1 − 2 ⋅ r6 ⋅ 
+
+ 1   ⋅ Tmi +1,6 +
kI
 
 kI
  
 2 ⋅ r6 ⋅ E6 ⋅ hre, 6 
 2 ⋅ r6 ⋅ E6 ⋅ he ,6 
 ⋅ TCÉU + 
 ⋅ T∞
+ 
kI
kI




(3.52)
Nodos superficiais internos
Nesta seção são apresentadas as equações algébricas que caracterizam os nodos
superficiais internos de cada parede da carroceria. As linhas tracejadas da Figura 3.13
representam os contornos dos seis planos que constituem as superfícies internas de cada
parede.
Figura 3.13 – Planos representativos das superfícies internas das paredes
Cap. 3 – Modelagem matemática
58
Para o nodo superficial interno da parede número 1, a energia que aflui ao mesmo
provém do interior da parede por transferência de calor por condução, e do ar ambiente
interno por transferência de calor por convecção.
O balanço de energia no nodo superficial interno da parede número 1 conduz à forma
explícita da equação em diferenças finitas, definida por:
 
h ⋅E 

h ⋅E
 
Tmi +,11 =  2 ⋅ r1 ⋅ i ,1 1  ⋅ Tmi −1 + 1 − 2 ⋅ r1 ⋅  i ,1 1 + 1  ⋅ Tmi ,1 +
kI 
 

 kI
 
+ (2 ⋅ r1 ) ⋅ T
i +1
m +1,1
(3.53)
onde:
r1
- critério de estabilidade para o nodo no instante i ;
hi ,1
- coeficiente de transferência de calor por convecção, W / m 2 ⋅ K ;
Tmi −1 - temperatura do ar interno, K ;
Tmi ,1 - temperatura superficial interna, parede 1, K ;
Tmi ++11,1 - temperatura superficial externa, parede 1, K;
obs.: para índice subscrito, i indica superfície interna e 1 indica parede;
para índice sobrescrito, i e i + 1 indicam o instante.
Para os nodos superficiais internos das paredes 2 a 6, a aplicação do balanço de energia
é semelhante àquele para superfície interna da parede número 1. A energia aflui ao nodo
interno provinda do interior da parede por transferência de calor por condução, e do ar
ambiente interno por transferência de calor por convecção.
Respeitados os índices pertinentes a cada parede, o balanço para o nodo superficial
interno da parede 2 é definido como:
 
h ⋅E 

h ⋅E
 
Tmi +, 21 =  2 ⋅ r2 ⋅ i , 2 2  ⋅ Tmi −1 + 1 − 2 ⋅ r2 ⋅  i , 2 2 + 1  ⋅ Tmi , 2 +
kI 
 

 kI
 
+ (2 ⋅ r2 ) ⋅ T
i +1
m +1, 2
Para o nodo superficial interno da parede 3, tem-se:
(3.54)
Cap. 3 – Modelagem matemática
 
h ⋅E 

h ⋅E
 
Tmi +,31 =  2 ⋅ r3 ⋅ i ,3 3  ⋅ Tmi −1 + 1 − 2 ⋅ r3 ⋅  i ,3 3 + 1  ⋅ Tmi ,3 +
kI 
 

 kI
 
+ (2 ⋅ r3 ) ⋅ T
i +1
m +1, 3
59
(3.55)
Para o nodo superficial interno da parede 4, tem-se:
 
h ⋅E 

 E ⋅h
 
Tmi +, 41 =  2 ⋅ r4 ⋅ i , 4 4  ⋅ Tmi −1 + 1 − 2 ⋅ r4 ⋅  4 i , 4 + 1  ⋅ Tmi , 4 +
kI 
 

 kI
 
+ (2 ⋅ r4 ) ⋅ Tmi ++11, 4
(3.56)
O nodo superficial interno da parede 5 apresenta:
 
h ⋅E 

h ⋅E
 
Tmi +,51 =  2 ⋅ r5 ⋅ i ,5 5  ⋅ Tmi −1 + 1 − 2 ⋅ r5 ⋅  i ,5 5 + 1  ⋅ Tmi ,5 +
kI 
 

 kI
 
+ (2 ⋅ r5 ) ⋅ Tmi ++11,5
(3.57)
O nodo da superficial interno da parede 6 está definido pela Equação 3.58.
h ⋅E 

h ⋅E
 
 
Tmi +,61 =  2 ⋅ r6 ⋅ i , 6 6  ⋅ Tmi −1 + 1 − 2 ⋅ r6 ⋅  i ,6 6 + 1  ⋅ Tmi , 6 +
kI 
 

 kI
 
+ (2 ⋅ r6 ) ⋅ Tmi ++11, 6
(3.58)
Ar ambiente interno
Para o balanço de energia são contabilizados os fluxos de energia provenientes das
superfícies internas das paredes da carroceria e da superfície da carga. A transferência
de calor por convecção ocorre através da superfície interna das paredes com o ar
ambiente interno, e deste, com a superfície da carga.
Na Figura 3.14, a linha contínua representa a carga inserida dentro da carroceria e a
linha tracejada representa a superfície interna das paredes da carroceria.
Cap. 3 – Modelagem matemática
60
Figura 3.14 – Disposição da carga no interior da carroceria
Fica estabelecido neste modelo que a transferência de calor por radiação entre as
superfícies internas das paredes não é contabilizada, devido aos reduzidos níveis de
temperatura envolvidos.
No modelo, a superfície inferior da carga fica apoiada sobre cavaletes dispostos sobre o
piso da carroceria. As demais superfícies da carga também ficam afastadas das
superfícies internas da carroceria. Com essa disposição, cria-se um espaço preenchido
com ar, o que contribui para uma melhora do isolamento da carga e favorece a
circulação do ar frio proveniente do sistema de refrigeração.
Aplicado o balanço de energia, a forma explícita da equação em diferenças finitas para
o ar ambiente interno é definida por:
i +1
m −1
T
 
 Ai ,1 ⋅ hi ,1 + Ai , 2 ⋅ hi , 2 + Ai ,3 ⋅ hi ,3 + Ai , 4 ⋅ hi , 4  
 
 ⋅ T i
= 1 − ri ⋅ Ci ⋅ 
+ A ⋅h + A ⋅h + h ⋅ A
  m−1

 
i
,
5
i
,
5
i
,
6
i
,
6
c
c

 
+ (ri ⋅ Ci ⋅ hi ,1 ⋅ Ai ,1 )⋅ Tmi +,11 + (ri ⋅ Ci ⋅ hi , 2 ⋅ Ai , 2 )⋅ Tmi +, 21 + (ri ⋅ Ci ⋅ hi ,3 ⋅ Ai ,3 )⋅ Tmi +,31
+ (ri ⋅ Ci ⋅ hi , 4 ⋅ Ai , 4 ) ⋅ Tmi +, 41 + (ri ⋅ Ci ⋅ hi ,5 ⋅ Ai ,5 )⋅ Tmi +,51 + (ri ⋅ Ci ⋅ hi , 6 ⋅ Ai , 6 )⋅ Tmi +, 61
+ (ri ⋅ Ci ⋅ hc ⋅ Ac ) ⋅ Tci − (ri ⋅ Ci ⋅ qe )
(3.59)
61
Cap. 3 – Modelagem matemática
onde:
ri
- critério de estabilidade para o nodo no instante i ;
Ai ,n - área superficial interna correspondente à parede n , m²;
Ac - área superficial da carga, m²;
hc - coeficiente de transferência de calor por convecção, W / m 2 ⋅ K ;
Tci - temperatura superficial da carga, K;
qe - fluxo de calor trocado com o evaporador, W/m²;
obs.: para índice subscrito, c indica superfície da carga;
para índice sobrescrito, i e i + 1 indicam o instante, n indica o número da
parede.
O coeficiente Ci é definido como:
Ci =
E pp
k i ⋅ Ai
(3.60)
onde:
E pp - espessura média da camada de ar ambiente interno, m;
Ai - área superficial interna total das paredes, m²;
Superfície da carga
O material que deve ser conservado sob refrigeração no interior da carroceria está
envolvido e troca calor com o ar ambiente interno. Ambos estão contidos pela superfície
interna das paredes e podem ser visualizados na Figura 3.14.
Para a determinação da temperatura da superfície da carga, deve-se avaliar o modelo de
transferência de calor em sólidos que melhor se aplica às condições existentes. Bejan
(1996) apresenta um estudo detalhado, evidenciando as faixas de aplicação do Método
da Capacitância Global e Método do Sólido Semi-infinito, para as condições de
62
Cap. 3 – Modelagem matemática
temperatura de superfície constante (caso 1), fluxo de calor constante na superfície
(caso 2) e superfície em contato com um fluido (caso 3).
O Método da Capacitância Global é utilizado quando o número de Biot (Bi) for inferior
a 0,1, conforme estudo apresentado por Incropera e Dewitt (1996). Este adimensional é
definido como o produto do coeficiente convectivo carga-ar e o comprimento
característico dividido pela condutividade térmica da carga, ou seja:
Bi =
hc ⋅ L
kc
(3.61)
onde:
hc - coeficiente de transferência de calor por convecção, W / m 2 ⋅ K ;
Neste caso, a equação para a discretização da temperatura da superfície da carga é
expressa na forma:
  2 ⋅ rc ⋅ hc ⋅ Ec
Tci +1 = 1 − 
kc
 
 i  2 ⋅ rc ⋅ hc ⋅ Ec
 ⋅ Tc + 
kc


 i +1
 ⋅ Tm −1

(3.62)
onde:
rc
- critério de estabilidade para o nodo no instante i ;
Ec - profundidade da carga, m;
k c - condutividade térmica da carga, W / m ⋅ K ;
obs.: para índice sobrescrito, i + 1 indica o instante. Para índice subscrito, c
indica superfície da carga, m − 1 indica ar interno;
Neste trabalho, foram implementadas apenas as equações referentes aos casos 2 e 3 do
Modelo para Sólido Semi-Infinito.
Para o caso 2, adotou-se a equação geral proposta por Bejan (1996) na forma:
Tc ( X , t ) − Ti = 2 ⋅
q´´  α ⋅ t 
⋅

k  π 
1/ 2


X 2  q´´
X
 −
⋅ exp  −
⋅ X ⋅ erfc 
1/ 2
 4 ⋅α ⋅ t  k
 2 ⋅ (α ⋅ t )




(3.63)
63
Cap. 3 – Modelagem matemática
que simplificada para a condição X = 0 se reduz a:
12
T0 (t ) = 2 ⋅
q"  α ⋅ t 
⋅

k  π 
(3.64)
+ Ti
onde:
q" - fluxo de calor, W/m²;
k - condutividade térmica do material, W / m ⋅ K ;
α - difusividade térmica do material, m²/s;
t - instante, s;
Ti - temperatura do inicial do corpo de prova, °C;
A forma discretizada é expressa como:
i +1
c
T
1


q"  α c ⋅ t  2 
i

= 2⋅ ⋅
 + Tc
 kc  π  


(3.65)
Para o caso 3, adotou-se a equação apresentada por Incropera e Dewitt (1996) na forma:
Tc ( X , t ) − Ti
 X
= erfc ⋅ 
(T∞ − Ti )
 2⋅ α ⋅t

h ⋅ α ⋅t
   hc ⋅ X hc2 ⋅ α ⋅ t  
X
 ⋅ erfc
 − exp
+
+ c


k
k
kc
   c
c
 
 2⋅ α ⋅t




(3.66)
Aplicando-se X = 0 na Equação 3.66, para obter a temperatura superficial da carga em
qualquer instante t , tem-se:

  h 2 ⋅ α ⋅ t  
 h ⋅ α ⋅t
 ⋅ erfc c
Tc (t ) =  3,9223 *10 −5 − exp c

kc

  k c  

(
)
 
  ⋅ (T∞ − Ti ) + Ti

 
(3.67)
64
Cap. 3 – Modelagem matemática
A forma discretizada da Equação 3.67 é:
Tc
i +1

 h ⋅ αc ⋅ t
  h 2 ⋅ α ⋅ t  
=  3,9223 *10 −5 − exp c c   ⋅ erfc c

kc
  

  k c

(
)
  i +1
  ⋅ Tm −1 − Tci

 
(
)
+ Tci
(3.68)
Resolução do sistema de equações
Para resolução do sistema de equações algébricas correspondentes à distribuição de
temperaturas nos nodos adotou-se o método de Gauss-Seidel. Este método iterativo
tem-se mostrado bastante eficiente para sistemas com grande número de equações
Maliska (1995).
No primeiro passo, são recalculadas as temperaturas em cada nodo a partir dos valores
previamente inicializados a partir das equações discretizadas. Este processo é repetido
até que o critério de convergência seja satisfeito para as temperaturas de todos os nodos
a cada novo instante i , ou de acordo com:
T j(i +1) − T j(i )
T j(i +1)
≤γ
(3.69)
onde:
T j(i +1) - temperatura no ponto j , no instante i + 1 ;
T j(i ) - temperatura no ponto j , no instante i ;
γ
- diferença relativa aceitável para temperatura, decimal;
Inicia-se, então, o cálculo das temperaturas para o instante seguinte, obedecendo os
mesmos procedimentos estabelecidos acima. O procedimento é mantido até que o tempo
total de simulação, definido como dado de entrada, seja atingido. Os fluxogramas das
Figuras 3.15 e 3.16 apresentam, respectivamente, o desenvolvimento destes dois passos.
65
Cap. 3 – Modelagem matemática
Início
Fornecer valores iniciais estimados
para as temperaturas
Entrar com valores para os diversos
coeficientes e propriedades dos materiais
e fluidos
1
Calcular o
temperaturas
novo
valor
Critério de
convergência é
atendido para todas
as temperaturas
encontradas?
S
das
N
Atualizar os
valores das
temperaturas
desconhecidas
com os novos
valores
encontrados
Temperatura
encontrada
1
2
Figura 3.15 – Fluxograma para o cálculo das temperaturas desconhecidas, passo 1
66
Cap. 3 – Modelagem matemática
2
Entrar com valores das temperaturas
para cada ponto, encontradas no
cálculo para o instante atual
Entrar com valores atualizados para os
diversos coeficientes e propriedades dos
materiais e fluídos, para temperaturas de
cada ponto, no instante atual
3
Calcular o valor de cada uma das
temperaturas desconhecidas, para
o instante atual
Critério de convergência
é atendido para todas as
temperaturas no instante
atual?
N
S
Passar ao
próximo
instante
2
N
Tempo total
decorrido ≥
instante final?
Atualizar os
valores das
temperaturas
desconhecidas
com os novos
valores
encontrados
S
3
S
Temperatura
final
para cada ponto
Fim
Figura 3.16 – Fluxograma para o cálculo das temperaturas desconhecidas, passo 2
Cap. 3 – Modelagem matemática
67
3.4- Decisão final
O procedimento final consiste em confrontar a capacidade que o conjunto motorequipamento de refrigeração selecionado possui, com a demanda de refrigeração
existente no interior da carroceria.
Para tal, ao final da simulação, o programa retorna um gráfico com duas curvas, onde
uma corresponde à capacidade de refrigeração oferecida pelo conjunto motorequipamento de refrigeração avaliados, e outra mostra o comportamento da carga
térmica instantânea que penetra no interior da carroceria, ao longo de um intervalo de
tempo predeterminado.
Portanto, o efeito de refrigeração oferecido pelo equipamento de refrigeração deve ser
necessariamente superior à carga de refrigeração obtida nos cálculos efetuados.
Caso tal premissa não seja atendida, os dados característicos de um novo modelo de
equipamento de refrigeração devem ser fornecidos ao programa para uma nova
avaliação, até que a premissa seja validada.
CAPÍTULO 4
VALIDAÇÃO DO MODELO E
DISCUSSÃO DOS RESULTADOS EXPERIMENTAIS
Neste capítulo são apresentados os procedimentos experimentais bem como os
resultados obtidos nos ensaios realizados, visando a validação do modelo matemático.
Foram realizados três ensaios experimentais distintos, onde dois cilindros, sendo um de
alumínio e outro de cobre e, também, uma peça de músculo bovino, foram inseridos em
uma estufa e ficaram expostos às condições ambientais internas controladas, durante o
intervalo de tempo determinado para cada ensaio.
Os ensaios foram realizados no laboratório de Transferência de Calor, nas dependências
da PUC Minas.
4.1- Materiais e métodos
Os ensaios foram realizados em uma estufa com 600 mm de aresta, mantida afastada do
piso e da parede, cujas superfícies interna e externa são construídas em chapa de aço.
As superfícies externas são pintadas na cor verde, conforme mostra a Figura 4.1a. O
aquecimento do ar interno é promovido por resistência elétrica, fixada sobre a face
interior de cada superfície interna, que está isolada termicamente por uma manta de lã
de rocha com espessura de 40 mm.
A parede superior possui uma abertura com tampa, mostrado em detalhe na Figura 4.1b,
para inserção dos corpos de prova.
68
69
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
Para controlar a potência elétrica fornecida à estufa foi utilizado um Varivolt, modelo
VM 215, com potência máxima de 1,5 kVA.
(a)
(b)
Figura 4.1 – Estufa – vista geral (a) e detalhe da abertura (b)
Para aferição da potência elétrica fornecida à estufa, foram utilizados um voltímetro
com capacidade máxima de 127 volts e resolução de 1Volt, além de um amperímetro
com capacidade máxima de 20 ampéres e resolução de 0,1 A, mostrados na Figura 4.2.
Figura 4.2 – Detalhe do conjunto Voltímetro-Amperímetro
Para obtenção das medidas de temperatura foram utilizados dois multímetros da marca
Polimed, modelo PM2700 e com resolução de 0,001 mV. O cronômetro da marca
Technos, modelo S08039, é digital e possui resolução de 0,01 s. As dimensões
geométricas foram determinadas com auxílio de uma trena da marca Starret, modelo
Tru-Lock, com resolução de 1 mm.
70
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
A Fig. 4.3 mostra o aparato experimental instalado e ligado, pronto para realização dos
ensaios.
Figura 4.3 – Aparato experimental utilizado nos ensaios
Corpos de prova
Os corpos de prova dos ensaios 1 e 2 são dimensional e construtivamente idênticos,
sendo maciços e diferindo somente em sua composição, conforme mostra a Fig. 4.4.
Ambos possuem 160 mm de altura, 50 mm de diâmetro e superfícies levemente
oxidadas.
(a)
(b)
Figura 4.4 – Corpos de prova - detalhe
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
71
Cada corpo de prova é posicionado verticalmente e montado a 200 mm da superfície
interna da parede interior de uma tampa idêntica à da estufa. Através de um orifício
feito longitudinalmente no corpo do cilindro, foi inserido um termopar tipo T, com saída
para conexão com o multímetro, conforme mostra a Figura 4.5.
Figura 4.5 – Detalhe construtivo do corpo de prova
Para o ensaio 3, foi utilizado como corpo de prova uma peça de músculo bovino, com
160 mm de altura, 50 mm de largura e 60 mm de profundidade, mostrada na Figura 4.6.
Figura 4.6 – Corpo de prova 3
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
72
Com formato aproximado de um paralelepípedo, foi fixada verticalmente à uma altura
aproximada de 200 mm da tampa. Buscou-se, dessa forma, reproduzir as condições
usadas nos ensaios 1 e 2.
Para efetuar a medida da temperatura na peça foram utilizados dois termopares tipo K.
Foi feito um orifício longitudinal no corpo da peça, onde o primeiro sensor foi
introduzido e posicionado à 80 mm da superfície. O segundo sensor foi fixado
superficialmente à peça. A Fig. 4.7 mostra o corpo de prova preparado para o ensaio 3.
Figura 4.7 – Corpo de prova posicionado para início do ensaio
Procedimento experimental
Com antecedência mínima de seis horas do início do ensaio, todo o aparato necessário
à experimentação foi montado. O varivolt e o voltímetro-amperímetro foram ligados à
rede elétrica e conectados aos cabos que fornecem eletricidade à estufa. Em seguida,
selecionou-se no varivolt a potência a ser fornecida ao conjunto.
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
73
Decorrido o intervalo de tempo necessário ao aquecimento e atingida a condição de
estabilidade da temperatura do ar interno, os termopares eram conectados aos
multímetros para leitura dos dados.
Com o cronômetro zerado, a tampa da estufa foi retirada e o corpo de prova inserido
rapidamente. O cronômetro era acionado no instante em que a tampa com o corpo de
prova fechava a estufa. A seguir, eram realizadas as leituras experimentais.
Tratamento dos dados
A partir dos dados obtidos nos multímetros, realizou-se o tratamento dos dados
experimentais a partir do polinômio recomendado em Controle e Instrumentação (1988,
1989), ou seja:
T = a 0 + a1 ⋅ E + a 2 ⋅ E 2 + a3 ⋅ E 3 + a 4 ⋅ E 4 + a5 ⋅ E 5 + a6 ⋅ E 6 + a7 ⋅ E 7 + a8 ⋅ E 8 (4.1)
onde:
T - temperatura, °C;
E - força eletromotriz, V;
a n - coeficientes do polinômio, mostrados na Tabela 4.1.
Tabela 4.1 – Coeficientes para termopar tipo T e tipo K
Coeficiente
a0
a1
a2
a3
a4
a5
a6
a7
a8
Valor
Tipo T
0,100860910
25727,94369
-767345,8295
78025595,81
-9247486589
6,97688E+11
-2,66192E+13
3,94078E+14
-
Tipo K
0,226584602
24152,10900
67233,42480
2210340,682
-860963914,9
4,83506E+10
-1,18452E+12
1,38690E+13
-6,33708E+13
4.2- Convecção Natural – Ensaio 1
Este ensaio consistiu basicamente em monitorar o aquecimento do cilindro de alumínio
no interior da estufa por 20 minutos. A potência elétrica responsável pelo aquecimento
do ar em seu interior foi fixada em 147 W.
74
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
As temperaturas iniciais do ar interno e do corpo de prova foram medidas instantes
antes do início do ensaio e durante sua realização. Os resultados obtidos são mostrados
nas Tabelas 4.2 e 4.3, respectivamente.
Tabela 4.2 – Temperaturas ambiente interna, externa e do corpo de prova iniciais.
Temperatura (°C)
Ar ambiente interno
109,2
Local
Ar ambiente externo
23,0
Corpo de prova
20,4
Tabela 4.3 – Evolução temporal das temperaturas experimentais para o Ensaio 1
Intervalo
(minutos)
0
1
2
3
4
5
6
7
8
9
10
12
14
16
18
20
C orpo de prova
Ar ambiente Interno
Valor leitura T em peratura Valor leitura T em peratura
(m V)
(°C )
(m V)
(°C )
-0,03
22,4
3,79
109,2
0,04
24,2
3,72
107,7
0,11
25,9
3,70
107,3
0,17
27,3
3,70
107,3
0,23
28,8
3,71
107,5
0,29
30,2
3,71
107,5
0,35
31,7
3,72
107,7
0,41
33,1
3,72
107,7
0,46
34,3
3,72
107,7
0,52
35,8
3,73
107,9
0,57
36,9
3,73
107,9
0,68
39,6
3,73
107,9
0,78
41,9
3,73
107,9
0,88
44,3
3,73
107,9
0,97
46,4
3,73
107,9
1,02
47,6
3,73
107,9
Validação do modelo matemático
Os cálculos experimentais indicaram que o número de Biot para o ensaio era da ordem
de 0,0044. Portanto, para validação do modelo matemático foi selecionada a rotina
correspondente ao Método da Capacitância Global. Sob estas condições, a distribuição
de temperatura espacial para a carga é uniforme, conforme Özisik (1990).
Como não foi possível medir as temperaturas internas e externas da estufa, foram
arbitrados como dados de entrada os primeiros valores recalculados na simulação
matemática e mostrados na Tabela 4.4.
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
75
Tabela 4.4 – Temperaturas iniciais das superfícies externas e internas
Temperatura
°C
°C
Parede
Superfície Externa Superfície Interna
Superior
30,4
113,0
Inferior
51,0
117,5
Lateral esquerda
30,4
113,0
Lateral direita
30,4
113,0
Frontal
30,4
113,0
Traseira
30,4
113,0
O tempo total de simulação foi definido em 20 minutos.
Considerações gerais
Para a seqüência de simulações realizadas, alguns parâmetros de entrada do programa
computacional foram ajustados de modo a reproduzir melhor as condições encontradas
no ensaio experimental, a saber:
-
Radiação solar: o valor da radiação incidente sobre a superfície
superior externa foi anulado.
-
Aquecimento da estufa: foi criado um termo adicional em que a
potência elétrica fornecida era uniformemente distribuída entre as seis
superfícies internas.
-
A temperatura ambiente externa e a temperatura do céu foram
consideradas constantes e iguais a 23,5 °C e 21,9 °C, respectivamente.
-
Ajuste da geometria: para o programa computacional desenvolvido a
forma geométrica adotada para a carga é a de um paralelepípedo, com a
maior dimensão na direção-y. Entretanto, no ensaio realizado em
laboratório o corpo de prova era cilíndrico, com a altura coincidente
com a direção –y. Portanto, nesta direção havia similaridade entre os
modelos. Considerou-se, então, o paralelepípedo de seção quadrada,
cujo lado foi determinado para se garantir volumes iguais entre os
corpos de prova. Os valores obtidos são mostrados na Tabela 4.5.
A Tabela 4.6 mostra as propriedades térmicas adotadas na simulação matemática.
76
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
Tabela 4.5 – Dimensões geométricas e propriedades térmicas do
alumínio adotadas na simulação matemática.
Propriedade
Dimensões Geométricas
Localização
Valor
Altura
150mm
Largura
40mm
Comprimento
40mm
Condutividade térmica
Alumínio
204 W/m.K*
Difusividade Térmica
Alumínio
8,42 x10-5 m2/s*
(*) Fonte: Özisik (1990)
Tabela 4.6 – Propriedades térmicas dos materiais utilizados na simulação
Propriedade
Emissividade da superfície
Localização
Valor
Externa – parede 1
0,95
Externa – paredes 2 a 6
0,95
Interna – paredes 1 a 6
0,88*
Condutividade térmica
Isolante
0,042 W/m.K
Difusividade Térmica
Isolante
3,50 x10-7 m2/s
(*) Para as paredes internas, foi considerada a emissividade para chapas de aço altamente oxidadas,
segundo Infrared –Thermography (2004).
Resultados Obtidos
Os resultados experimentais e simulados para a temperatura do corpo de prova são
mostrados na Tabela 4.7 e no gráfico da Figura 4.8.
77
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
Tabela 4.7 – Estudo comparativo entre os resultados experimentais
e simulados para o Ensaio 1
Temperatura (ºC)
T em po
(minutos)
0
2
3
4
5
6
7
8
9
10
12
14
16
18
20
T em peratura do corpo de prova
Valor experimental
Sim ulação
(°C )
(°C )
22,4
22,4
25,9
23,9
27,3
25,4
28,8
26,9
30,2
28,4
31,7
29,8
33,1
31,2
34,3
32,6
35,8
33,9
36,9
35,3
39,6
37,9
41,9
40,4
44,3
42,8
46,4
45,2
47,6
46,3
D esvio
Absoluto
(°C )
0,0
2,0
1,9
1,9
1,8
1,9
1,9
1,7
1,9
1,6
1,7
1,5
1,5
1,2
1,3
50,0
45,0
40,0
35,0
30,0
25,0
20,0
0
5
10
15
20
25
Tempo (minutos)
Experimental
Simulado
Figura 4.8 – Avaliação comparativa entre as temperaturas experimentais e
simuladas do corpo de prova – Ensaio 1
Tais resultados evidenciam que os desvios absolutos médio e máximo são iguais a 1,7 e
2,0ºC, respectivamente. Estas diferenças foram consideradas aceitáveis para validação
do modelo proposto, visto que os desvios esperados em medições com termopares tipo
T, para temperaturas inferiores a 400ºC, é da ordem de 3ºC e de 1ºC para calibrações
segundo a norma DIN 43710 e ANSI MC 96.1 – 1975, respectivamente, Controle e
Instrumentação (1988, 1989).
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
78
Deve-se destacar, ainda, que o modelo matemático desenvolvido não incluiu a troca
radiante de calor entre o corpo de prova e as paredes internas, devido aos níveis
reduzidos das temperaturas envolvidas. Entretanto, para as condições em que o Ensaio 1
foi realizado tal termo poderia ser incluído para melhorar os resultados obtidos na
simulação matemática.
4.3- Convecção Natural - Ensaio 2
O segundo ensaio é similar ao anterior, tendo sido utilizado como corpo de prova um
cilindro maciço de cobre. Para garantir o nível de aquecimento da estufa, a potência
elétrica foi mantida igual a 147 W. O tempo total de ensaio foi fixado em 21 minutos.
As temperaturas iniciais do ar interno e do corpo de prova foram medidas instantes
antes do início do ensaio e durante sua realização. Os resultados obtidos são mostrados
nas Tabelas 4.8 e 4.9, respectivamente.
Validação do modelo matemático
O número de Biot calculado para o ensaio 2 foi de 0,0023, justificando-se a utilização
do Método da Capacitância Global. Da mesma forma, todas as considerações discutidas
para o ensaio 1 foram aplicadas a este novo teste.
As temperaturas adotadas para inicialização do método numérico são mostradas na
Tabela 4.10.
Os valores adotados para simulação matemática são os mesmos apresentados nas
Tabelas 4.5 e 4.6, exceto as propriedades térmicas típicas para o cobre. A condutividade
térmica do cobre puro é de 398 W/m.K e sua difusividade térmica igual a 1,16x10-4m2/s,
Incropera e Dewitt (1998).
Tabela 4.8 – Temperaturas iniciais do ar ambiente interno, externo e do corpo de prova
Temperatura (°C)
Ar ambiente interno
110,3
Local
Ar ambiente externo
24,0
Corpo de prova
23,7
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
79
Tabela 4.9 – Evolução temporal das temperaturas experimentais para o Ensaio 2
Intervalo
(minutos)
0
1
2
3
4
5
6
7
8
9
10
11
12
13
15
17
19
21
C orpo de prova
Ar ambiente Interno
Valor leitura T em peratura Valor leitura T em peratura
(m V)
(°C )
(m V)
(°C )
-0,02
23,7
3,80
110,3
0,03
24,9
3,73
108,8
0,08
26,1
3,72
108,6
0,13
27,3
3,72
108,6
0,17
28,3
3,72
108,6
0,22
29,5
3,72
108,6
0,27
30,7
3,72
108,6
0,32
31,9
3,73
108,8
0,36
32,9
3,72
108,6
0,41
34,1
3,73
108,8
0,46
35,3
3,74
109,0
0,50
36,2
3,74
109,0
0,55
37,4
3,74
109,0
0,59
38,4
3,74
109,0
0,68
40,5
3,73
108,8
0,77
42,6
3,73
108,8
0,85
44,5
3,73
108,8
0,93
46,4
3,73
108,8
Tabela 4.10– Temperaturas iniciais das superfícies
T em peratura
°C
°C
Parede
Superfície Externa Superfície Interna
Superior
30,5
114,0
Inferior
51,0
118,5
Lateral esqu erda
30,5
114,0
Lateral direita
30,5
114,0
Frontal
30,5
114,0
T raseira
30,5
114,0
Resultados Experimentais
Os resultados experimentais e simulados para a temperatura do corpo de prova no
ensaio 2 são mostrados na Tabela 4.11 e no gráfico da Figura 4.9.
Para esse ensaio, os desvios absoluto médio e máximo são iguais a 2,3 e 3,5ºC,
respectivamente. Os maiores desvios observados podem ser explicados pelo
desconhecimento do nível de pureza do cobre utilizado no corpo de prova, a ausência do
termo radiante nos balanços de energia, termopares não calibrados, dentre outros
fatores.
80
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
Tabela 4.11 – Resultados experimentais comparados aos resultados da simulação 2
T em peratura do corpo de prova
Valor experimental
Sim ulação
(°C )
(°C )
23,7
23,7
26,1
24,8
27,3
25,8
28,3
26,9
29,5
27,9
30,7
28,9
31,9
30,0
32,9
31,0
34,1
32,0
35,3
32,9
36,2
33,9
37,4
34,8
38,4
35,8
40,5
37,6
42,6
39,4
44,5
41,2
46,4
42,9
T em po
(minutos)
0
2
3
4
5
6
7
8
9
10
11
12
13
15
17
19
21
D esvio
Absoluto
(°C )
0,0
1,3
1,5
1,4
1,6
1,8
1,9
1,9
2,1
2,4
2,3
2,6
2,6
2,9
3,2
3,3
3,5
Temperatura (ºC)
5 0 ,0
4 5 ,0
4 0 ,0
3 5 ,0
3 0 ,0
2 5 ,0
2 0 ,0
0
5
10
15
20
25
T e m p o ( m in )
E x p e rim e n t a l
S im u la d o
Figura 4.9 – Avaliação comparativa entre as temperaturas experimentais e
simuladas do corpo de prova – Ensaio 2
4.4- Convecção Natural + Radiação – Ensaio 3
Neste ensaio, uma peça composta por músculo bovino e com formato aproximado ao de
um paralelepípedo, foi submetida a aquecimento monitorado dentro da estufa durante 72
minutos. O aquecimento do ar interno é proporcionado pelo fornecimento à estufa de
57W de potência elétrica.
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
81
As temperaturas para o ar ambiente interno e externo e, também, as temperaturas
superficial e interna do corpo de prova são mostradas nas Tabelas.4.12 e 4.13.
Tabela 4.12 – Temperaturas iniciais para o ar ambiente e para o corpo de prova
Tempo
(minutos)
0
Temperatura do ar ambiente
Interno
Externo
(°C)
(°C)
52,1
23,5
Temperatura do corpo de prova
Superfície
Interna
(°C)
(°C)
20,1
19,8
Tabela 4.13 – Evolução temporal das temperaturas do corpo de prova experimentais –
Ensaio 3
T em po
(minutos)
0
2
4
6
9
12
15
19
23
27
32
37
42
47
52
57
62
67
72
Superfície externa
T em peratura
Valor leitura
(mV)
(°C )
-0,12
20,1
-0,09
20,8
-0,06
21,5
-0,04
22,0
0,00
23,0
0,04
24,0
0,07
24,7
0,12
25,9
0,16
26,9
0,20
27,8
0,24
28,8
0,29
30,0
0,33
31,0
0,37
32,0
0,41
33,0
0,45
33,9
0,48
34,7
0,52
35,6
0,58
36,4
Interior
T em peratura
Valor leitura
(m V)
(°C )
-0,13
19,8
-0,13
19,8
-0,13
19,8
-0,13
19,8
-0,12
20,1
-0,11
20,3
-0,08
21,0
-0,04
22,0
-0,02
22,5
0,01
23,2
0,05
24,2
0,12
25,9
0,17
27,1
0,20
27,8
0,26
29,3
0,30
30,3
0,35
31,5
0,39
32,5
0,42
33,3
Validação do modelo matemático
Os cálculos experimentais indicaram que o número de Biot para o ensaio 3 era da ordem
de 2. Assim, adotou-se o caso 2 do Método do Sólido Semi- infinito, com o objetivo de
se incluir o termo de troca de calor radiante.
O tempo total de simulação foi definido em 57 minutos.
82
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
Novamente, foram arbitrados como dados de entrada os primeiros valores recalculados
na simulação matemática e mostrados na Tabela 4.14.
Tabela 4.14 – Temperaturas iniciais das superfícies
T em peratura
°C
°C
Parede
Superfície Externa Superfície Interna
Superior
25,3
53,4
Inferior
34,5
55,7
Lateral esqu erda
25,3
53,4
Lateral direita
25,3
53,4
Frontal
25,3
53,4
T raseira
25,3
53,4
Considerações gerais
Para a seqüência de simulações realizadas, as mesmas considerações discutidas para os
ensaios 1 e 2 foram mantidas. Os valores obtidos para os ajustes das dimensões
geométricas da carga são mostrados na Tabela 4.15.
Tabela 4.15 – Dimensões geométricas e propriedades térmicas da carne
adotadas na simulação matemática.
Propriedade
Localização
Dimensões Geométricas Altura
Valor
160mm
Largura
50mm
Comprimento
60mm
Condutividade térmica
Carne bovina
0,48 W/m.K*
Difusividade Térmica
Carne bovina
1,32x10-7 m2/s*
(*)Fonte: Singh (1998)
Para aplicação do caso 2 do Método do Sólido Semi-infinito, o fluxo de calor total na
superfície, convectivo e radiante, deve ser constante durante todo o processo. Neste
estudo, foram realizados cálculos preliminares em que foi verificada a variação do fluxo
total de energia durante o período a ser avaliado. Assim, para emprego do método, o
83
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
fluxo transiente foi decomposto como o somatório de fluxos constantes durante
pequenos intervalos de tempo, com boa aproximação.
Resultados Experimentais e Simulados
Os resultados obtidos são mostrados na Tabela 4.16 e no gráfico da Figura 4.10
Tabela 4.16 - Resultados experimentais comparados aos resultados da simulação 3
T em po
(minutos)
0
2
4
6
9
12
15
19
23
27
32
37
42
47
52
57
T em peratura do corpo de prova
Valor experim ental
Simulação
(°C )
(°C )
20,1
20,1
20,8
23,0
21,5
23,8
22,0
25,1
23,0
25,9
24,0
26,6
24,7
27,2
25,9
28,0
26,9
28,6
27,8
28,9
28,8
29,6
30,0
30,0
31,0
30,5
32,0
30,8
33,0
31,3
33,9
31,6
D esvios
Absoluto
(°C )
0,0
2,2
2,3
3,1
2,9
2,6
2,5
2,1
1,7
1,1
0,8
0,0
-0,5
-1,2
-1,7
-2,3
O desvio absoluto médio encontrado foi de 1ºC. Entretanto, analisando-se os dados
obtidos, verifica-se que o desvio absoluto sofre um rápido crescimento durante os 6
primeiros minutos. A partir desse instante, o desvio absoluto é reduzido, voltando a
crescer após os 37 minutos. Nesta primeira fase, os valores simulados são sempre
superiores ao medidos.
Atribui-se essa resposta lenta da temperatura simulada na superfície da carne a fatores
não previstos no modelo matemático, como alteração das propriedades físicas e térmicas
da carne com o aumento de sua temperatura superficial e às perdas de calor por
evaporação decorrentes de seu elevado teor de umidade inicial.
84
Temperatura (ºC)
Cap. 4 – Validação do modelo e discussão dos resultados experimentais
40,0
30,0
20,0
10,0
0,0
0
20
40
60
Tempo (minutos)
Tsup - exp
Tsup - sim
Figura 4.10 – Comparação dos resultados experimentais e simulados da
Temperatura superficial da carne – Ensaio 3.
CAPÍTULO 5
SIMULAÇÃO MATEMÁTICA - ANÁLISE DE
RESULTADOS
Neste capítulo é apresentado um estudo de caso sobre o transporte e armazenamento de
carne bovina, que deve ser mantida resfriada no interior de um baú refrigerado, com
objetivo de demonstrar a aplicação da metodologia a uma situação prática.
Inicialmente, são definidas as condições iniciais, dimensões e parâmetros de projeto do
baú frigorífico da referida simulação. Tais condições buscam refletir um ensaio em
campo, realizado na cidade de Campo Belo/ MG. Infelizmente, condições restritivas
impostas pelo próprio armazenamento comprometeram a fase de medições,
principalmente devido à impossibilidade de efetuá-las em todos os pontos internos do
baú.
5.1- Estudo de caso – Resfriamento de carne
O baú frigorífico estava montado em um caminhão, devidamente carregado e
estacionado ao ar livre.
Materiais e métodos
O ensaio de campo foi realizado em um pátio, onde estava estacionado o furgão
frigorífico selecionado para o ensaio junto a dois outros caminhões, conforme mostrado
na Figura 5.1. Nessa figura, também pode ser vista a porta traseira de acesso ao seu
interior.
O furgão foi previamente condicionado, sendo depois carregado com a carne préresfriada. As medições foram realizadas enquanto o mesmo aguardava a ordem de
viagem.
85
Cap. 5 – Simulação matemática – Análise de resultados
86
Figura 5.1 – Vista dos caminhões estacionados
O carregamento de carne é composto de peças traseiras, dianteiras e costelas bovinas,
dependuradas no teto, conforme mostrado na Figura 5.2.
Figura 5.2 – Disposição da carga
Este furgão era equipado com um equipamento de refrigeração convencional, em que o
compressor principal funciona ligado diretamente ao motor do veículo.
Durante o período em que o motor do caminhão ficava desligado, a refrigeração era
proporcionada por um compressor auxiliar, acoplado ao motor elétrico que permanecia
conectado à rede de energia elétrica, conforme detalhe na Figura 5.3.
Cap. 5 – Simulação matemática – Análise de resultados
87
Figura 5.3 – Caixa de tomadas para fornecimento de energia elétrica
O acionamento desse motor é controlado por um termostato posicionado no interior da
carroceria, que obedece aos limites máximo e mínimo definidos como aceitáveis para a
temperatura do ar interno.
As medições das temperaturas internas foram efetuadas após o desligamento automático
do sistema auxiliar de refrigeração. Em seguida, foram realizadas as medições das
condições externas (temperaturas das superfícies e ambiente e velocidade do vento).
Todas as temperaturas das superfícies interna e externa e da própria carga resfriada
foram medidas com o termômetro digital infra-vermelho da marca Omega, modelo OS
651, com resolução de 1°C. Para tal, eram ajustadas as emissividades relativas a cada
superfície e material empregado.
As medidas das temperaturas dos ambientes interno e externo e da velocidade do vento
foram obtidas à partir de um termo-anemômetro digital da marca Omega, modelo HHF
710. Este instrumento opera com resoluções de 0,1°C e 0,01 m/s, respectivamente para
medidas de temperatura e velocidade.
As dimensões da carroceria foram obtidas antecipadamente, antes do carregamento ser
efetuado. Foi utilizada para isso uma trena marca Starret, modelo Tru-Lock, com
resolução de 1mm.
Cap. 5 – Simulação matemática – Análise de resultados
88
Quantificação da energia a ser disponibilizada
O primeiro passo é a caracterização do motor e respectivos dados operacionais. Neste
estudo, foi utilizado um motor diesel 6L com ciclo de quatro tempos e turboalimentado.
Os dados de desempenho desse motor foram obtidos através de um dinamômetro e são
apresentados por Horuz (1999).
Nas Figuras 5.4 e 5.5 apresenta-se o comportamento da vazão mássica para diversos
regimes de carga e o perfil de temperaturas dos gases de exaustão, ambos em função da
potência útil do motor.
Figura 5.4 – Vazão dos gases de exaustão x potência útil do motor.
Adaptado de Horuz (1999)
Figura 5.5 – Temperatura dos gases de exaustão x potência útil do motor
Adaptado de Horuz (1999)
89
Cap. 5 – Simulação matemática – Análise de resultados
Para essa simulação, os valores típicos que foram adotados para estes parâmetros
operacionais são mostrados na Tabela 5.1 e são válidos para um regime de carga
equivalente a 200 Nm de torque e 30 kW de potência de saída, conforme Horuz (1999).
Tabela 5.1 – Parâmetros operacionais do motor
Parâmetro
Vazão dos gases de exaustão
(kg/s)
Temperatura dos gases de exaustão (°C)
Valor
0,1
300
Para o motor diesel, Heywood (1988) recomenda utilizar a razão de equivalência
combustível/ar da ordem de 0,8, que caracteriza uma faixa estequiométrica pobre.
Seleção do equipamento de refrigeração por absorção
Adotou-se, nessa simulação, o ciclo por absorção por simples efeito, que utiliza
água/amônia como fluido de trabalho. Conforme a Tabela 2.1, para o funcionamento
deste ciclo, torna-se necessário a disponibilização de uma fonte térmica com
temperatura entre 120-150 °C. O coeficiente de performance é considerado igual a 0,5.
A efetividade do trocador de calor para recuperação da energia contida nos gases de
exaustão foi arbitrada como 0,5.
Quantificação dos fluxos de energia para o interior do baú
Para o desenvolvimento do balanço de energia para o baú frigorífico devem ser
definidas inicialmente a localização e condições ambientais externas.
Foi adotado um dia de céu claro no mês de janeiro para a cidade de Campo Belo. O
caminhão encontrava-se estacionado em local plano e aberto, com a ocorrência limitada
de ventos.
Para esta localidade, os parâmetros de simulação como localização geográfica e
condições ambientais externas, são mostrados nas Tabela 5.2 e 5.3, respectivamente.
Tabela 5.2 – Parâmetros geográficos da localidade
Parâmetro
Latitude
Longitude
Altitude (m)
Valor
21° 14' S
45° 00'
920
90
Cap. 5 – Simulação matemática – Análise de resultados
Tabela 5.3 – Dados de entrada adotados da simulação
Parâmetro
Total de dias do mês
Dia do ano
Horário de início da simulação
Umidade relativa do ar*
Temperatura máxima para o mês*
Temperatura mínima para o mês*
Velocidade do vento
Velocidade do veículo
Valor
31 dias
17
13:30h
81,3%
27,8 °C
17,7 °C
0,8 m/s
0 m/s
Os dados assinalados com (*) foram extraídos das Normais Climatológicas, DNMET
(1992). A velocidade local do vento média foi estimada em 0,8 m/s a partir de uma
série de medidas experimentais realizadas no local.
Caracterizações superficiais e do interior das paredes
O baú frigorífico possui o formato de um paralelepípedo, posicionado horizontalmente
sobre o chassi do caminhão, cujos parâmetros dimensionais e construtivos são
apresentados na Tabela 5.4.
Tabela 5.4 – Parâmetros dimensionais e construtivos do baú
Parâmetro
Dimensões externas
Comprimento
Largura
Altura
Espessura do isolamento
Parede 1
Paredes de 2 a 6
Valor
7,50 m
2,60 m
2,40 m
0,10 m
0,05 m
O isolamento térmico é de espuma rígida de poliuretano, com condutividade térmica da
ordem de 0,023 W / m ⋅ K , conforme manual da empresa Kingspan Insulation. Sua
difusividade térmica foi calculada, sendo o valor encontrado igual a 3,5 x 10-7 m2/s.
Externamente, as paredes laterais, frontal e traseira do baú são revestidas por chapas
frisadas de alumínio anodizado. O teto é revestido por chapas lisas de alumínio
anodizado e o piso é constituído por um painel compensado de madeira.
Os parâmetros relativos à transferência de calor por radiação para os materiais que
compõem as superfícies externas das paredes do baú estão mostrados na Tabela 5.5.
91
Cap. 5 – Simulação matemática – Análise de resultados
Tabela 5.5 – Propriedades óticas dos materiais e coberturas das superfícies externas
Parâmetro
Emissividade da superfície externa da parede 1
Absortividade da superfície externa da parede 1
Emissividade da superfície externa das paredes 3 a 6
Valor
0,44
0,44
0,44
As superfícies externas e internas das paredes, essas últimas revestidas com placas em
P.R.F.V, foram consideradas isotérmicas.
A Tabela 5.6 mostra as temperaturas superficiais das paredes medidas, utilizadas como
valores iniciais na simulação. A temperatura da superfície interna da parede 3 foi
estimada, devido à impossibilidade de alcançar o fundo da carroceria.
Tabela 5.6 – Temperaturas iniciais das superfícies externas e internas
Parede/ Nº
Superior/1
Inferior/2
Traseira/6
Frontal/5
Lateral esquerda/3
Lateral direita/4
Temperatura
°C
°C
Superfície
Superfície
54,0
14,0
32,0
11,0
51,0
11,0
25,0
11,0
28,0
11,0
26,0
11,0
Condições ambientais internas
O baú carrega uma carga de carne resfriada com formato de um paralelepípedo, que está
posicionado horizontalmente em seu interior. As dimensões desta carga foram
estimadas e correspondem a um bloco com as dimensões apresentadas na Tabela 5.7.
Tabela 5.7 – Parâmetros dimensionais da carga
Parâmetro
Comprimento
Largura
Altura
Valor
7,2 m
1,3 m
1,3 m
O ventilador e o equipamento de refrigeração permaneceram desligados. Com isso a
condição de movimentação do ar no interior do baú foi considerada estacionária, sendo
a quantidade de energia trocada com o evaporador nula.
A temperatura inicial do ar ambiente interno era igual a 10°C e a da carga era de 9 °C.
92
Cap. 5 – Simulação matemática – Análise de resultados
As propriedades térmicas para a carga foram obtidas em Singh (1998), e são mostradas
na Tabela 5.8.
Tabela 5.8 – Parâmetros para a carga
Carga
Temperatura média da superfície da carga
Condutividade térmica
Difusidade térmica
Valor
9° C
0,48 W/mK
1,32 x 10E-7
O programa computacional contabiliza apenas a radiação que incide sobre o teto do baú.
Os parâmetros referentes ao calor gerado por essa radiação constam da Tabela 5.9.
Os valores para os coeficientes de Bennett são referentes ao mês de janeiro, conforme
Pereira (2002). O valor da insolação diária em média mensal é válido para a cidade de
Lavras (MG) e foram obtidos a partir das Normais Climatológicas, DNMET (1992).
Tabela 5.9 – Parâmetros relacionados a radiação solar absorvida
Parâmetro
Coeficientes empíricos da equação de Bennett
a
b
c
Insolação diária em média mensal (horas)
Valor
0,225
0,4812
0,0007
187,9
Considerações
Nesta simulação matemática, a temperatura da superfície da carga foi obtida adotandose o modelo do sólido semi-infinito com superfície em contato com um fluído (caso 3,
conforme apresentado no Capítulo 4).
Embora a carne estivesse acondicionada em peças individuais, no modelo matemático
adotou-se a forma de um paralelepípedo. Para determinação de seu volume, tomou-se a
massa total de carne embarcada (13.125 kg) que dividida por sua densidade (1.086
kg/m3) fornecia o volume equivalente, da ordem aproximada de 12 m3.
Resultados Obtidos
O tempo de simulação foi fixado em 25 minutos, cujos resultados para as temperaturas
das superfícies externas das paredes são apresentados na Tabela 5.10 e no gráfico da
Figura 5.6.
93
Cap. 5 – Simulação matemática – Análise de resultados
Tabela 5.10 – Temperaturas superficiais externas das paredes
Tempo
(s)
0
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
Text 1
(°C)
54,0
55,0
55,1
55,1
55,1
55,1
55,1
55,0
55,0
54,9
54,8
54,8
54,7
54,6
54,5
54,4
54,2
54,1
54,0
53,9
53,7
53,6
53,5
53,3
53,2
Text 2
(°C)
32,0
28,2
27,9
27,6
27,3
27,1
26,9
26,8
26,7
26,6
26,5
26,5
26,4
26,4
26,3
26,3
26,2
26,2
26,2
26,2
26,1
26,1
26,1
26,1
26,1
Text 3
(°C)
51,0
26,4
25,6
25,1
24,8
24,5
24,4
24,3
24,2
24,2
24,1
24,1
24,1
24,1
24,1
24,1
24,1
24,1
24,1
24,1
24,1
24,0
24,0
24,0
24,0
Text 4
(°C)
25,0
24,1
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
24,0
Text 5
(°C)
28,0
24,0
23,8
23,6
23,5
23,5
23,4
23,4
23,4
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
Text 6
(°C)
26,0
23,7
23,6
23,5
23,4
23,4
23,4
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
23,3
Figura 5.6 – Gráfico das temperaturas superficiais externas das paredes
Cap. 5 – Simulação matemática – Análise de resultados
94
Analisando-se o comportamento das temperaturas das superfícies externas, verifica-se
uma pequena diminuição ao longo do tempo da ordem de 2ºC para as superfícies do
piso e da porta traseira. Para as outras superfícies, a redução do valor das temperaturas é
ainda menor. No teto, o fator preponderante é a incidência da radiação solar, que
justifica os elevados níveis de temperatura encontrados.
Deve-se salientar, ainda, que o modelo não avaliou corretamente a temperatura da
superfície traseira (Text3), conforme pode ser constatado nas duas primeiras linhas da
Tabela 5.10, porque este prevê a incidência de radiação solar somente sobre o teto.
Durante os ensaios, foi constatada a incidência de radiação solar nessa superfície de
forma significativa.
A evolução temporal das temperaturas da superfície da carga, do ar interno e das
superfícies internas das paredes são apresentados na Tabela 5.11 e no gráfico da Figura
5.7.
Tabela 5.11 – Temperaturas superficiais internas das paredes, do ar interno e da carga
Tempo Tint 1 Tint 2 Tint 3 Tint 4 Tint 5 Tint 6 Tar int Tcarga
(s)
(°C) (°C) (°C) (°C) (°C) (°C)
(°C)
(°C)
0
14,0 11,0 11,0 11,0 11,0 11,0
10,0
9,0
2
14,7 12,5 13,7 12,1 12,2 12,1
11,4
9,0
3
14,8 12,7 13,7 12,3 12,4 12,3
11,3
9,0
4
14,9 12,9 13,8 12,5 12,5 12,4
11,4
9,1
5
15,0 13,0 13,8 12,6 12,6 12,6
11,5
9,1
6
15,1 13,2 13,8 12,8 12,8 12,7
11,6
9,1
7
15,3 13,3 13,8 12,9 12,9 12,8
11,7
9,1
8
15,4 13,4 13,9 13,0 13,0 13,0
11,7
9,1
9
15,5 13,5 13,9 13,1 13,1 13,1
11,8
9,1
10
15,6 13,6 13,9 13,2 13,2 13,2
11,9
9,1
11
15,7 13,7 14,0 13,3 13,3 13,3
11,9
9,1
12
15,7 13,7 14,0 13,3 13,3 13,3
11,9
9,1
13
15,8 13,7 14,0 13,4 13,4 13,4
12,0
9,1
14
15,9 13,8 14,0 13,5 13,5 13,4
12,0
9,1
15
16,0 13,9 14,1 13,6 13,6 13,5
12,1
9,2
16
16,1 13,9 14,1 13,7 13,6 13,6
12,1
9,2
17
16,2 14,0 14,1 13,8 13,7 13,7
12,2
9,2
18
16,2 14,1 14,2 13,8 13,8 13,7
12,2
9,2
19
16,3 14,1 14,2 13,9 13,8 13,8
12,3
9,2
20
16,4 14,2 14,2 14,0 13,9 13,9
12,3
9,2
21
16,5 14,2 14,3 14,0 14,0 13,9
12,3
9,2
22
16,6 14,3 14,3 14,1 14,0 14,0
12,4
9,2
23
16,6 14,3 14,3 14,2 14,1 14,0
12,4
9,2
24
16,7 14,3 14,4 14,2 14,1 14,1
12,4
9,2
25
16,8 14,4 14,4 14,3 14,2 14,2
12,5
9,2
Cap. 5 – Simulação matemática – Análise de resultados
95
Figura 5.7 - Temperaturas superficiais internas das paredes, do ar interno e da carga
A análise das temperaturas das superfícies internas mostra um comportamento bastante
similar. O aumento de suas temperaturas verificado durante o período de simulação,
decorrente da transferência de calor por condução através das paredes, variou na faixa
de 0,8 a 2,1 ºC.
Conforme era esperado, a temperatura superficial interna do teto (Tint1) é superior à das
outras superfícies a despeito da maior espessura do isolamento térmico adotado porque
a temperatura externa era bastante elevada.
Para a temperatura do ar interno, constata-se uma elevação de 2,0 °C acima de seus
valores iniciais, evidenciando as trocas térmicas com as superfícies internas das paredes,
enquanto que a temperatura da carga sofreu uma pequena variação durante os 25
minutos.
É importante salientar que, na prática, o equipamento de refrigeração opera conforme a
indicação dos sensores da temperatura do ar interno do baú. De acordo com a faixa de
temperaturas pré-selecionada nos termostatos, o intervalo entre o desligamento e o
religamento do equipamento de refrigeração pode ser menor que o intervalo de tempo
definido nesta simulação.
Cap. 5 – Simulação matemática – Análise de resultados
96
O comportamento das temperaturas do céu e do ar ambiente externo, assim como a
energia absorvida pela superfície externa do teto decorrente da radiação solar incidente
( q s ) foram mostrados na Figura 5.8.
Figura 5.8 – Temperaturas do céu e do ar ambiente externo e energia solar
absorvida pela superfície do teto do baú frigorífico
Análise Final
Na simulação matemática desenvolvida, deve-se avaliar se a capacidade de refrigeração
requerida para atender as condições iniciais pré-definidas pode ser suprida pelo sistema
de refrigeração por absorção.
Na Figura 5.9, a curva que representa o comportamento da carga térmica no interior da
carroceria é representada pela linha contínua. A linha tracejada representa a capacidade
de refrigeração do sistema selecionado. Observa-se que o efeito de refrigeração
desejado atende durante todo o período de simulação à demanda requerida.
Cap. 5 – Simulação matemática – Análise de resultados
97
Figura 5.9 – Capacidade de refrigeração disponibilizada x demanda por refrigeração
CAPÍTULO 6
CONCLUSÕES E CONSIDERAÇÕES FINAIS
Neste trabalho, foi desenvolvida uma metodologia para análise da viabilidade de
aproveitamento da energia contida nos gases de exaustão de motores diesel, através de
um equipamento de refrigeração por absorção, para condicionamento de produtos
perecíveis em caminhões baú.
A metodologia proposta se baseia na avaliação comparativa entre a carga térmica
requerida, dependente do tipo de produto, temperatura de armazenamento e condições
ambientais típicas encontradas durante o transporte, e a energia disponível nos gases de
exaustão usada como aporte energético no ciclo de refrigeração por absorção.
Essa metodologia, discutida no Capítulo 3, foi implementada em um programa
computacional que permite simular cargas diversas sob diferentes condições de
operação. Para tal, foram desenvolvidas rotinas para aplicação dos modelos
matemáticos em sólidos semi-infinitos e método da capacitância global.
Ensaios realizados em laboratório permitiram a validação das diferentes rotinas
desenvolvidas. Os resultados foram considerados satisfatórios com desvios absolutos
máximos entre os valores das temperaturas experimentais e simuladas inferior a 2,5ºC.
Tais resultados conduzem às seguintes conclusões:
•
o aproveitamento dos gases de exaustão de um motor, para geração de efeito de
refrigeração através de um sistema de absorção é viável, para o caso estudado;
•
o programa computacional, desenvolvido neste trabalho, calcula a parcela útil da
energia térmica contida nos gases de exaustão de um motor diesel, que pode ser
utilizada para gerar efeito de refrigeração em um equipamento de absorção;
98
Cap. 6 – Conclusões e considerações finais
•
99
o programa computacional, com algumas restrições e ajustes prevê com razoável
precisão a evolução temporal das temperaturas do ar interno e das superfícies
internas e externas das paredes de uma câmara com superfícies internas isotérmicas;
•
o modelo matemático desenvolvido apresenta desvios mais significativos entre
valores experimentais e simulados, quando estas envolvem aquecimento de material
orgânico com alto teor de umidade;
Feitas estas colocações, sugere-se para trabalhos futuros a inclusão na modelagem
matemática do termo relativo à troca radiante entre as superfícies internas das paredes,
alteração da rotina de caracterização da carga, dando maior flexibilidade à definição dos
parâmetros geométricos e às especificidades de transporte, como o uso de embalagens
individuais.
Em termos práticos, sugere-se o desenvolvimento e montagem de uma bancada
experimental de um sistema de refrigeração por absorção para trabalhar em conjunto
com um motor Diesel em teste. Dessa forma, o efeito de refrigeração produzido no
evaporador do ciclo pode ser avaliado em função dos diferentes regimes de operação do
motor.
REFERÊNCIAS BIBLIOGRÁFICAS
Bibliografia Consultada
ASHRAE, 1993, Fundamentals. American Society of Heating, Refrigeration and Airconditioning Engineers, Inc., Atlanta.
ASHRAE, 1995, Technical Data Bulletin. Desiccant and Absorption Cooling. v.11,
n.2. American Society of Heating, Refrigeration and Air-conditioning Engineers, Inc.,
Atlanta.
ASM INTERNATIONAL, 1997, Handbook Committee. Materials Handbook. USA,
Metals Parks.
Anônimo ,“Princípios e aplicações práticas de termopares”. Revista Controle &
Instrumentação, nov., 1988, pág. 32-38, dez., 1988, pág. 17-21, jan 1989, pág. 20-25.
ANTAS, Luiz Mendes, 1980, Dicionário de termos técnicos: inglês - português, 6 ed.,
São Paulo, Traço editora.
BEJAN, Adrian, 1995, Convection heat transfer, 2 ed., USA, John Wiley & Sons, inc.
BEJAN, Adrian, 1996, Transferência de calor, 1 ed., Tradução de Euryale de Jesus
Zerbini e Ricardo Santilli Ekman Simões, São Paulo, Edgard Blücher Ltda.
BURGHARDT, M. David, 1982, Engineering thermodynamics with applications, 2 ed.,
New York, Harper & Row.
BURMEISTER, Louis C., 1983, Convective heat transfer, 1 ed., USA, John Wiley &
Sons Inc.
CARVALHO, José Guilherme de, 1991, “Alternativas para uso do gás natural em
sistemas de refrigeração”. Revista Abrava, março/abril, pág. 30 – 33.
CHAPMAN, Alan J., 1984, Heat transfer, 4 ed., New York, Collier Macmillan Canada,
Inc.
100
101
CHINNAPPA, J. C. V., 1992, “Principles of Absorption Systems Machines” in Sayigh,
A. A. M. e McVeigh, J.C., editors, Solar air conditioning and refrigeration, 1 ed.,
Oxford, Pergamon Press, pág. 13 – 65.
COOPER, Jeffery, 1998, Introduction to partial differential equations with MATLAB,
USA, Birkhäuser.
CORTEZ, Luís Augosti B., MÜHLE, Ingo Norberto, e SILVA, Andrés da , 1994,
“Refrigeração por absorção com o par água–amônia e seu potencial no caso brasileiro”.
Revista Abrava, jan./fev., pág. 33-38.
DNMET, 1992, Normais climatológicas : 1961-1990, Brasília.
DUFFIE, John A., e BECKMAN, William A., 1991, Solar engineering of thermal
processes, 2 ed., Canada, John Wiley & Sons, Inc.
FELAMINGO, José Carlos, 2001, “O uso de refrigeração por absorção com
recuperação de calor”. Revista do frio, fev e março, pág. 48 – 52 e pág. 44 – 48.
FIGLIOLA, Richard S. e BEASLEY, Donald E., 1991, Theory and design for
mechanical measurements, 1 ed., Singapore, John Wiley & Sons, Inc.
GANDER, Walter e HREBÍCEK, Jirí, 1995, Solving problems in cientific computing
using Maple and MATLAB, 2 ed., Berlin, Germany, Springer-Verlag Hidelberg New
York.
GARCIA, Cláudio, 1997, Modelagem e simulação de processos industriais e de
sistemas eletromecânicos, 1 ed., São Paulo, Editora da Universidade de São Paulo.
GILES, Ranald V., Mecânica dos fluidos e hidráulica, Tradução de Sérgio dos Santos
Borde, São Paulo, McGraw-Hill do Brasil.
GIVONI, Baruch, 1994, Passive and low energy cooling of buildings, USA, John
Wiley & Sons, Inc.
HARTNETT, James P. e IRVINE Jr, Thomas F., 1990, Advances in heat transfer, vol.
20, California, Academic Press, Inc.
102
HARTNETT, James. P. e ROHSENOW, Warren. M., 1998, Handbook of heat transfer,
3 ed., USA, McGraw-Hill.
HEROLD, Keith. E., RADERMACHER, Reinhard e KLEIN, Sanford. A., 1996,
Fundamentos de unidades resfriadoras e bombas de calor por absorção. Tradução de
Roberto A. Peixoto do livro Absorption Chillers and Heat Pumps, CRC Press, Inc.,
Seminário de Refrigeração por Absorção, São Paulo, 2001.
HEYWOOD, John B., 1988, Internal combustion engine fundamentals, 1 ed., New
York, McGraw – Hill.
HORUZ, Ilhami., 1999, “Vapor absorption refrigeration in road transport vehicles”,
Journal of energy engineering, August/1999, pág. 48 – 58.
HOLMAN, Jack Philip, 1983, Transferência de calor, Tradução de Luiz Fernando
Milanez, São Paulo, McGraw-Hill do Brasil.
HUANG, Francis F., 1982, Engineering thermodynamics: fundamentals and
applications, 2 ed., Macmillan.
INCROPERA, Frank P. e DEWITT, David P., 1998, Fundamentos de Transferência de
Calor e de Massa, 4 ed., Tradução de Sérgio Stamile Soares, Rio de Janeiro, LTC
editora.
JABARDO, José M. Saiz. “Amônia em sistemas frigoríficos”. Revista Abrava, jan/fev,
1994, p. 17 – 32.
JABARDO, José M. Saiz. “Refrigerantes”. Tecnologia da refrigeração, n. 7, 2001, p.
22 –29.
JOAQUIM, Ana Paula. “Fluido refrigerante: a aplicação ideal para cada tipo de
instalação”. Tecnologia da refrigeração, n. 7, 2001, pág. 12 – 21.
KAVIANY, Massoud, 1994, Principles of convective heat transfer, 1 ed., New York,
Springer-Verlag New York, Inc.
KERN, Donald Q., 1980, Processos de Transmissão de calor, 1 ed., Tradução de Adir
M. Luiz, Rio de Janeiro, Guanabara Dois S.A.
103
KHARAB, Abdelwahab e GUENTHER, Ronald B., 2002, An introduction to numerical
methods: a MATLAB approach, Flórida, CRC Press.
KNUDSEN, James G. e KATZ, Donald L., 1958, Fluid dynamics and heat transfer, 1
ed., New York, McGraw-Hill Book company, inc.
KREITH, Frank, 1977, Princípios da transmissão de calor, 3 ed., Tradução por Eitaro
Yamane e outros, São Paulo, Edgard Blücher Ltda.
LEWIS, Roland W. e MORGAN, Kenneth, 1981, Numerical methods in heat transfer,
1 ed., USA, John Wiley & Sons Ltd.
LI, Kam W., 1996, Applied thermodynamics: availability method and energy
conversion, New York, Taylor & Francis.
LINDBORG, Anders. “Amônia: referência e avaliação de risco quando utilizada como
refrigerantes em sistemas de refrigeração”. Revista Abrava, jan./fev., 1996, pág. 32–40.
MALISKA, Clovis R., 1995, Transferência de calor e mecânica dos fluidos
computacional: fundamentos e coordenadas generalizadas, Rio de Janeiro, LTC
editora.
MATLAB, The Language of Technical Computing, Version 6.5, 2002, The
MathWorks, Inc.
MCADAMS, William H., 1954, Heat transmission, 3 ed., New York, McGraw-Hill
Book Company Inc.
MCQUISTON, Faye C. e PARKER, Jerald D., 1994, Heating, ventilating, and air
conditioning: analysis and design, 4 ed., John Wiley & Sons, Inc.,
MIKHAILOV, Mikhail Dimitrov, e ÖZIŞIK, M. Necati, 1993, Unified analysis and
solutions of heat and mass diffusion, Canada, General Publishing Company, Ltda.
ÖZIŞIK, M. Necati, 1990, Transferência de Calor: um texto básico, 1 ed., Tradução de
Luiz de Oliveira. Rio de Janeiro, Editora Guanabara.
ÖZIŞIK, M. Necati, 1994, Finite Difference Methods in Heat Transfer, 1 ed., Flórida,
CRC Press, Inc.
104
PAO, Yen-Ching, Engineering analysis: interactive methods and programs with
FORTRAN, QuickBASIC, MATLAB, and Mathematica, New York, CRC Press.
PATANKAR, Suhas V., 1980, Numerical heat transfer and fluid flow, Series in
computation methods in mechanics and thermal sciences, 1 ed., USA, Hemisphere
Publishing Corporation.
PEREIRA, Elizabeth Marques Duarte, 1997, Siscos, Cuba.
PEREIRA, Elizabeth Marques Duarte, 2001, Energia solar térmica: instalações solares
de pequeno porte. Belo Horizonte, PUC Minas Virtual.
PEREIRA, Elizabeth Marques Duarte et al., 2002, Energia solar aplicada: instalações
solares de pequeno porte. Belo Horizonte, PUC Minas Virtual.
RADHAKRISHNAN, Sudhaharini, 1997, Measurement of thermal properties of
seafood, Dissertação de mestrado, Biological Systems Engineering Department,
Virginia Polytechnic Institute and State University, Blacksburg, Virginia.
REINICK, Amélia Rubíolo de, 1994, “Avaliação de parâmetros para estimar os tempos
de congelamento e descongelamento”, Revista Abrava, jul./ago., pág. 40-44.
RUGGIERO, Márcia A. Gomes, 1996, Cálculo numérico: aspectos teóricos e
computacionais. São Paulo, Makron Books.
SANTOS, Jurandir Augusto dos, 2000, Óleo Diesel: produtos Petrobrás, 5 ed., Rio de
Janeiro, Petrobrás.
SBRAVATI, Alan, 2001, “Projeto de sistema para condicionamento de ar por absorção
de água-brometo de lítio”, Revista do Frio, pág. 41-52.
SCHILICHTING, Hermann, 1979, Boundary-layer theory, 4 ed., New York, McGrawHill Company.
SHAMES, Irving Herman, 1973, Mecânica dos fluidos, v.1, Princípios básicos, 1 ed.,
Tradução de Mauro O. C. Amorelli, São Paulo, Edgard Blücher Ltda.
SIEGEL, Robert e HOWELL, John R., 1992, Thermal radiation heat transfer, 3 ed.,
USA, Publishing Office.
105
SILVA, Andrés da, 2001, “Avaliação energética e exergética de um sistema de
refrigeração por absorção água-amônia”, Revista do Frio, fev., pág. 36-46.
SINGH, R. Paul e HELDMAN, Dennis R., 1993, Introducción a la ingeniería de los
alimentos, 2 ed., Espanha, Zaragoza, Editorial Acribia.
SMITH, Charles C., LÖF, George, e JONES, Randy, 1994, Measurement and analysis
of evaporation from na inactive outdoor swimming pool, Solar Energy, v. 53, n. 1, pág.
3-7, USA, Elsevier Science Ltd.
SONNTAG, Richard E., BORGNAKKE, Claus, PARK, Kyoung Kuhn e PARK, Young
Moo, 1996, Catt2, Computer-Aided Thermodynamics Tables 2, Version 1.0a, John
Wiley & Sons, Inc.
SONNTAG, Richard E. e BORGNAKKE, Claus, 2003, Introdução à termodinâmica
para engenharia, Tradução de Luiz Machado, Geraldo Augusto Campolina França e
Ricardo Nicolau Nassar Koury, Rio de Janeiro, LTC Editora.
THOMAS, Lindon C., 1985, Fundamentos de transferência de calor, Tradução de
Priscila Aya Shimizu Günther e Cláudio Augusto Oller do Nascimento, Rio de Janeiro,
Prentice-Hall do Brasil Ltda.
WARK, Kenneth., 1977, Thermodynamics, 3 ed., London, MacGraw Hill.
WHITE, Frank M., 1994, Fluid mechanic, 3 ed., New York, McGraw-Hill, Inc.
WYLEN, Gordon J. Van e SONNTAG, Richard E., 1976, Fundamentos da
Termodinâmica clássica, 2 ed., Tradução de Eitaro Yamane e outros, São Paulo, Edgard
Blücher Ltda.
WYLEN, Gordon J. Van, SONNTAG, Richard E. e BORGNAKKE, Claus, 2000,
Fundamentos da termodinâmica, Tradução da 5ª edição americana por Euryale de Jesus
Zerbini, São Paulo, Edgard Blücher Ltda.
YU, Wen-Shing, LIN, Hsiao-Tsung, HWANG, Tsung-Yuan, 1991, “Conjugate heat
transfer of conduction and forced convection along wedges and a rotating cone”, Int. J.
Heat Mass Transfer, v. 34, n. 10, pág. 2.497 – 2.507.
106
Bibliografia Complementar
DIEHL, Peter, HAUBNER, Frank, KLOPSTEIN, Stefan, and KOCH, Franz, 2001,
“Exhaust Heat Recovery System for Modern Cars”. SAE Technical paper series.
Disponível em: < http://www.sae.org/technical/papers/2001-01-1020 >
GORDON, Sanford e MCBRIDE, Bonnie J., 1994, “Computer Program for Calculation
of Complex Chemical Equilibrium Compositions and Applications”. NASA Reference
Publication 1311
Disponível em: < http://gltrs.grc.nasa.gov/reports/1994/RP-1311.pdf >
LOWENSTEIN, Andrew, SLAYZAK, Steve. J., RYAN, Joseph P., PESARAN, Ahmad
A., 1998, “Advanced Commercial Liquid-Desiccant Technology Development Study”,
National Renewable Energy Laboratory, NREL/TP-550-24688.
Disponível em: < http://www.nrel.gov/docs/fy99osti/24688.pdf >
MCBRIDE, Bonnie J., ZEHE, Michael J. and GORDON, Sanford, 2002, “NASA Glenn
Coefficients for Calculating Thermodynamic Properties of Individual Species”. NASA /
TP – 2002-211556.
Disponível em: < http://gltrs.grc.nasa.gov/reports/2002/TP-2002-211556.pdf >
RAFFERTY, Kevin D., 1998, “Absortion Refrigeration”, Geo-Heat Center Quarterly
Bulletin, v. 19, n. 1, c. 13, Pág. 299 - 306.
Disponível em: < http://geoheat.oit.edu/pdf/tp51.pdf >
SLAYZAK, Steven. J., PESARAN, Ahmad. A., e HANCOCK, C. E., 1996,
“Experimental Evaluation of Commercial Desiccant Dehumidifier Wheels”, Nikanpour,
D. and Kosatte, S. Eds. Ab-sorption 96. Procedings international absorption heat pump
conference ’96, Quebec, Canada, September 17-20.
Disponível em: < http://www.nrel.gov/desiccantcool.wheels.html >
SRIKHIRIN, Pongsid, APHORNRATANA, Satha e CHUNGPAIBULPATANA,
Supachart, 2001, “A Review of Absorption Refrigeration Technologies”, in Renewable
and Sustainable Energy Reviews, Vol. 5, n. 4, p. 343-372.
Disponível em: < http://www.elsevier.com/locate/rser >
Electro Optical Indrustries, Inc., “Material Emissivity Properties”.
Disponível em: < http://www.electro-optical.com/bb_rad/emissivity/matlemisivty.asp >
107
Infrared Services, Inc., “Emissivity Values for Commom Materials”.
Disponível em: < http://www.infrared-thermography.com/material.htm >
Microtherm Thermal Insulation, “Thermal Conductivity Table”
Disp. em: < http://www.cla5ponline.org/download/energg_testing/2000/117/taiwan_paper.pdf >
Mikron Instrument Company, “Table of Emissivity of Various Surfaces”,
Disponível em: < http://www.mikron.com.br >
Square One Research PTY LTD, “Materials Properties of Solar Absorptance and
Emmittance”
Disponível em: < http://www.squ1.com >
Anexo A.1 – Listagem do programa
108
%========================================================================//
%// CTERMICA
//
%// Programa:
Autor: Valbert Garcia Assumpçao
//
%//
//
%//
Daniel Teixeira Gervasio
//
%//
Vinicius Meireles Ciriaco
//
%//
//
%//
//
%// Mestrado em Engenharia Mecanica - Pontificia Universidade Catolica de Minas Gerais - Brasil //
//
//
%//
//
%//=======================================================================//
clc
clear all;clear
%======================= Dados dos gases de exaustao do motor =====================
clc
massa=input(' Informe a massa de ar dos gases (kg/s): ');
T=input(' Informe a temperatura dos gases de saida do motor : ');
T=T+273.15;
fi=input(' Informe a razao de equivalencia (combustivel/ar): ');
R=8314.3; % Constante universal dos gases (Joule/kMol*Kelvin)
nCO2=12;
nH2O=13;
nO2=(18.5*(1-fi))/fi;
nN2=69.8/fi;
nT=nCO2+nH2O+nO2+nN2;
xCO2=nCO2/nT;
xH2O=nH2O/nT;
xO2=nO2/nT;
xN2=nN2/nT;
%------------------ Massa molecular ----------------------------------mCO2=44.01;
mH2O=18.015;
mO2=31.999;
mN2=28.013;
%-------------- Calculo de Cp para cada especie ------------------------cpCO2=(R/mCO2)*(((4.943650540*10^4)*(T^-2))+((-6.264116010*10^2)*(T^-1))+5.301725240+...
((2.503813816*10^-3)*T)+((-2.127308728*10^-7)*(T^2))+((-7.689988780*10^-10)*(T^3))...
+((2.849677801*10^-13)*(T^4)));
cpH2O=(R/mH2O)*(((-3.947960830*10^4)*(T^-2))+((5.755731020*10^2)...*(T^
1))+(9.317826530*10^-1)+((7.222712860*10^-3)*T)+((-7.342557370*10^6)*(T^2))+((4.955043490*10^-9)*(T^3))+((-1.336933246*10^-12)*(T^4)));
cpO2=(R/mO2)*(((-3.425563420*10^4)*(T^-2))+((4.847000970*10^2)*(T^-1))+1.119010961...
+((4.293889240*10^-3)*T)+((-6.836300520*10^-7)*(T^2))+((-2.023372700*10^-9)*(T^3))...
+((1.039040018*10^-12)*(T^4)));
cpN2=(R/mN2)*(((2.210371497*10^4)*(T^-2))+((-3.818461820*10^2)*(T^-1))+6.082738360...
+((-8.530914410*10^-3)*T)+((1.384646189*10^-5)*(T^2))+((-9.625793620*10^-9)*(T^3))...
+((2.519705809*10^-12)*(T^4)));
cpTot=((xCO2*cpCO2)+(xH2O*cpH2O)+(xO2*cpO2)+(xN2*cpN2)) % (Joule/kg*Kelvin);
Anexo A.1 – Listagem do programa
qEX=massa*cpTot*(T-423.15);
disp('
')
disp('
')
109
% (Watts)
disp([' Energia disponivel nos gases de exaustao: ', num2str(qEX), ' W' ]);
disp('
')
disp('
')
efet=input(' Informe a efetividade do trocador de calor ');
disp('
')
qG=qEX*efet;
% (Watts)
COP=input(' Informe o coeficiente de performance do equipamento de refrigeraçao ');
disp('
')
qEV=qG*COP;
disp('
% (Watts)
')
disp([' Efeito de refrigeraçao no evaporador: ', num2str(qEV), ' W' ]);
disp('
')
disp('
')
disp('
')
%================================= Latitude ==================================
disp(' LATITUDE')
graus=input(' Digite os graus: ');
if (graus>90)|(graus<-90)
while (graus>90)|(graus<-90)
errordlg('Dado incorreto!','LATITUDE, Graus')
disp(' Dado incorreto!')
graus=input(' Digite os graus novamente: ');
disp(' ')
end
end
minutos=input(' Digite os minutos: ');
if (minutos>60)|(minutos<-60)
while (minutos>60)|(minutos<-60)
errordlg('Dado incorreto!','LATITUDE, Minutos')
disp(' Dado incorreto!')
minutos=input(' Digite os minutos novamente: ');
disp(' ')
end
end
disp(' ')
latitude=graus+minutos/60;
%================================ Longitude =================================
disp(' LONGITUDE')
graus=input(' Digite os graus: ');
Anexo A.1 – Listagem do programa
110
if (graus>360)|(graus<0)
while (graus>360)|(graus<0)
errordlg('Dado incorreto!','LONGITUDE, Graus')
disp(' Dado incorreto!')
graus=input(' Digite os graus novamente: ');
disp(' ')
end
end
minutos=input(' Digite os minutos: ');
if (minutos>60)|(minutos<0)
while (minutos>60)|(minutos<0)
errordlg('Dado incorreto!','LONGITUDE, Graus')
disp(' Dado incorreto!')
minutos=input(' Digite os minutos novamente: ');
disp(' ')
end
end
disp(' ')
longitude=graus+minutos/60;
%================================ Altitude ==================================
disp(' ALTITUDE')
altitude=input(' Informe a altitude em metros: ');
altitude=altitude/1000;
disp(' ')
%============================= Umidade Relativa ===============================
disp(' UMIDADE RELATIVA DO AR')
umidader=input(' Digite o valor da umidade relativa do ar em(%): ');
if (umidader>100)|(umidader<1)
while (umidader>100)|(umidader<1)
errordlg('Dado incorreto!','UMIDADE RELATIVA DO AR')
disp(' Dado incorreto!');
umidader=input(' Digite o dado novamente: ');
disp(' ')
end
end
disp(' ')
umidader=umidader/100;
%==================================Mes e dia=================================
disp(' TEMPORAIS')
mes=input(' Informe a quantidade de dias no mes (28 - 29 - 30 - 31): ');
if (mes~=28)|(mes~=29)|(mes~=30)|(mes~=31)
while (mes~=28)|(mes~=29)|(mes~=30)|(mes~=31)
disp(' Dado incorreto!')
mes=input(' Digite novamente a quantidade de dias no mes (28 - 29 - 30 - 31): ');
disp(' ')
end
end
disp(' ')
dia=input(' Informe o dia (0 ate 365): ');
if (dia>366)|(dia<1)
while (dia>366)|(dia<1)
errordlg('Dado incorreto!','Dia do ano')
disp(' Dado incorreto!')
dia=input(' Digite novamente o dia (0 ate 365): ');
disp(' ')
end
Anexo A.1 – Listagem do programa
111
end
disp(' ')
hora=input(' Informe a hora (HH): ');
if (hora>24)|(hora<0)
while (hora>24)|(hora<0)
errordlg('Dado incorreto!','Hora')
disp(' Dado incorreto!')
hora=input(' Digite novamente a hora (0 ate 24): ');
disp(' ')
end
end
disp(' ')
%================================== Temperaturas =============================
disp(' TEMPERATURAS')
tmax=input(' Informe a temperatura maxima: ');
tmin=input(' Informe a temperatura minima: ');
if tmin>=tmax
while tmin>=tmax
errordlg('Temperatura minima deve ser menor que a maxima!','TEMPERATURAS')
disp(' Temperatura minima deve ser menor que a maxima!')
tmin=input(' Informe a temperatura minima: ');
disp(' ')
end
end
%================================= Hora Solar ================================
hs=hora;
deltat=tmax-tmin;
tambiente=tmax-(deltat/2)+(deltat/2)*cos((15*((hs)-14)*pi)/180)
pressaosat=7*10^(-8)*tambiente^3+7*10^(-8)*tambiente^2+tambiente*6*10^(-5)+0.0005;
pressaov=pressaosat*umidader;
%================================ Calculo do TDP =============================
pressaovv=pressaov*1000;
tdp=0.0065*pressaovv^5 - 0.1595*pressaovv^4 + 1.5601*pressaovv^3 - 7.9701*pressaovv^2+25.739...
*pressaovv - 12.192;
%========================== Calculo da Temperatura do Ceu ========================
tambiente=tambiente+273.15;
tceu=tambiente*(0.8+(tdp/250))^0.25;
tempceu=tceu-273.15;
disp([' Temperatura do ceu : ', num2str(tempceu)])
%================================= Carroceria ================================
disp(' ')
disp(' CALCULO DA CARGA TERMICA PARA A CARROCERIA')
disp(' ')
vvento=input(' Informe a velocidade do vento em m/s : ');
vveiculo=input(' Informe a velocidade do veiculo em km/h: ');
vveiculo=vveiculo/3.6;
disp([' Velocidade do veiculo em m/s: ', num2str(vveiculo)])
%========================== Definindo a Velocidade Relativa =======================
if (vveiculo==0)
vrelat=vvento;
else if (vvento>vveiculo)
Anexo A.1 – Listagem do programa
112
vrelat=vvento;
else vrelat=vveiculo;
end
end
%=========================================================================
disp(' ')
disp(' CALCULO DOS COEFICIENTES DE TRANSFERENCIA DE CALOR ')
disp('
POR CONVECÇAO ATRAVES DAS PAREDES')
disp(' ')
%======================== Calculo das dimensoes e Area Total =======================
disp(' DIMENSOES DA CARROCERIA')
Le=input(' Informe o comprimento externo da carroceria: ');
We=input(' Informe a largura externa da carroceria: ');
He=input(' Informe a altura externa da carroceria: ');
disp(' ')
Ae12=We*Le;
disp(' Soma das Areas externas 1 e 2')
disp([' ',num2str(Ae12)])
disp(' ')
Ae34=We*He;
disp(' Soma das Areas externas 3 e 4')
disp([' ',num2str(Ae34)])
disp(' ')
Ae56=Le*He;
disp(' Soma das Areas externas 5 e 6')
disp([' ',num2str(Ae56)])
disp(' ')
Aetotal=(Ae12+Ae34+Ae56)*2;
disp(' Area Total da Carroceria')
disp([' ',num2str(Aetotal)])
disp(' ')
%======================== Dimensoes Internas da Carroceria ========================
for i=1:6
EspI(i)=input([' Informe a espessura para a parede ', num2str(i),' (mm): ']);
EspI(i)=EspI(i)./1000;
disp([' ',num2str(EspI(i)),' metros'])
end
Wi=We-(EspI(5)+EspI(6));
Hi=He-(EspI(1)+EspI(2));
Li=Le-(EspI(3)+EspI(4));
fprintf('\n')
%=========================== Calculo das Areas Internas ==========================
Ai12=Wi*Li;
disp(' Soma das Areas internas 1 e 2')
disp([' ',num2str(Ai12)])
disp(' ')
Ai34=Wi*Hi;
disp(' Soma das Areas internas 3 e 4')
disp([' ',num2str(Ai34)])
disp(' ')
Ai56=Li*Hi;
disp(' Soma das Areas internas 5 e 6')
disp([' ',num2str(Ai56)])
disp(' ')
Ai=(Ai12+Ai34+Ai56)*2;
disp(' Area Interna')
disp([' ',num2str(Ai)])
Anexo A.1 – Listagem do programa
113
disp(' ')
%============= Calculo do coeficiente transferencia de calor por radiaçao externo =============
emissividadeext1=input(' Informe o valor da emissividade da superficie externa 1: ');
emissividadeext=input(' Informe o valor da emissividade para demais superficies externas: ');
absortext1=input(' Informe a absortividade da superficie externa 1: ');
tal=5.6697*10^(-8);
disp(' ')
%=========================== PAREDES EXTERNAS ============================
for n=1:6
disp(' ')
disp('Parede 1: teto ')
disp('Parede 2: piso ')
disp('Parede 3: lateral direita ')
disp('Parede 4: lateral esquerda ')
disp('Parede 5: traseira ')
disp('Parede 6: frente ')
fprintf('\n')
Text(n)=input([' Informe a temperatura superficial externa da parede ',num2str(n),': ']);
Text(n)=Text(n)+273.15;
%=========================================================================
disp(' ');
if (vrelat==0)
x3=0;
y3=1;
if (n==1)|(n==2)
z3=1;
s3=0;
else
z3=0;
s3=1;
end
else
x3=1;
y3=0;
if (n==1)|(n==2)|(n==5)|(n==6)
z3=1;
s3=0;
else
z3=0;
s3=1;
end
end
L(n)=Le*(x3)*(z3)+(We/2)*(x3)*(s3)+((Le*We)/((2*We)+(2*Le)))*(y3)*(z3)+He*(y3)*(s3);
tfluext(n)=(Text(n)+tambiente)/2;
betae(n)=1/tfluext(n);
%=============================== Calculo do Ro ===============================
ro(n)=351.92*tfluext(n)^(-1.0017);
%========================= Calculo do Cp, K e Alfa Externo ========================
cpexterno(n)=(0.103409*10)+(-0.28488708*10^(-3)*tfluext(n))+(0.7816818*10^(-6)*tfluext(n)^2)...
+(-0.4970786*10^(-9)*tfluext(n)^3)+(0.1077024*10^(-12)*tfluext(n)^4);
kexterno(n)=(-2.276501*10^(-3))+(1.2598485*10^(-4)*tfluext(n))+(-1.4815235*10^(7)*tfluext(n)^2)...
+(1.73550646*10^(-10)*tfluext(n)^3)+(-1.066657*10^(-13)*tfluext(n)^4)...
Anexo A.1 – Listagem do programa
114
+(2.47663035*10^(-17)*tfluext(n)^5);
alfaexterno(n)=(kexterno(n)/(ro(n)*cpexterno(n)))/1000;
%==================== Calculo da Viscosidade Dinamica e Cinematica ===================
viscosidadedin(n)=(-9.8601*10^(-1))+(9.080125*10^(-2)*tfluext(n))+(-1.17635575*10^(-4)...
*tfluext(n)^2)+(1.2349703*10^(-7)*tfluext(n)^3)+(-5.7971299*10^(-11)*tfluext(n)^4);
viscosidadecin(n)=viscosidadedin(n)/(ro(n)*1000000);
%============================= Calculo de Prandtl ==============================
prandtl(n)=viscosidadecin(n)/alfaexterno(n);
%============================= Calculo de Reynolds =============================
reynext(n)=vrelat*L(n)/viscosidadecin(n);
%=========================================================================
if (Text(n)>tambiente)
tx(n)=Text(n)-tambiente;
else
tx(n)=tambiente-Text(n);
end
%============================= Calculo de Rayleigh =============================
ra(n)=((9.8*betae(n)*tx(n)*(L(n)^3))/(viscosidadecin(n)*alfaexterno(n)));
%=========================== Teste da velocidade relativa ==========================
if (vrelat==0)
x1=0;
if (n==1)
y1=0;
z1=0;
s1=1;
else if (n==2)
y1=0;
z1=1;
s1=0;
end
if (n==3)|(n==4)|(n==5)|(n==6)
y1=1;
z1=0;
s1=0;
end
end
disp(' ')
%=============================== Calculo de Nusselt ============================
if (vrelat==0) & (Text(n)>tambiente)
x1=0;
y1=0;
if (n==1)
z1=1;
s1=0;
else
z1=0;
s1=1;
end
Anexo A.1 – Listagem do programa
115
if (n==2)
z1=0;
s1=1;
else
z1=1;
s1=0;
end
if (n==3)|(n==4)|(n==5)|(n==6)
y1=1;
z1=0;
s1=0;
end
end
else
x1=1;y1=0;z1=0;s1=0;
end
nuext(n)=((0.037*(reynext(n)^(4/5)))*(prandtl(n)^(1/3)))*x1...
+(((0.825+((0.387*(ra(n)^(1/6)))/((1+((0.492/prandtl(n))^(9/16)))^(8/27))))^2))*y1...
+(0.27*(ra(n)^(1/4)))*s1+(0.15*(ra(n)^(1/3)))*z1;
%============================ Impressao de dados ===============================
disp([' Superficie externa, face: ', num2str(n)])
disp(' ')
disp('
RO
Cp
K
Alfa
V ')
disp([ro(n);cpexterno(n);kexterno(n);alfaexterno(n);viscosidadecin(n)]')
fprintf('\n')
disp('
Prandtl Reynolds
Rayleigh Nusselt ')
disp([prandtl(n);reynext(n);ra(n);nuext(n)]')
disp(' ')
%=========== Coeficiente externo de transferencia de calor por convecçao e radiacao ==========
he(n)=(nuext(n)*kexterno(n))/L(n);
disp([' Coeficiente externo de convecçao vale ',num2str(he(n))]);
end
%=========================================================================
hre1=(Ae12)*emissividadeext1*tal*(tceu^2+Text(1)^2)*(tceu+Text(1));
hre3=(Ae34)*emissividadeext*tal*(tceu^2+Text(3)^2)*(tceu+Text(3));
hre4=(Ae34)*emissividadeext*tal*(tceu^2+Text(4)^2)*(tceu+Text(4));
hre5=(Ae56)*emissividadeext*tal*(tceu^2+Text(5)^2)*(tceu+Text(5));
hre6=(Ae56)*emissividadeext*tal*(tceu^2+Text(6)^2)*(tceu+Text(6));
%=========================================================================
fprintf('\n')
disp(' Coeficiente externo de transferencia de calor por radiaçao para a face 1');
hrext1=hre1/Ae12;
disp(hrext1)
%================================FIM EXTERNO==============================
%===============================INICIO INTERNO ============================
disp(' ')
disp(' CALCULO DOS COEFICIENTES DE TRANSFERENCIA DE CALOR POR CONVECÇAO')
disp(' PARA A SUPERFICIE INTERNA DA PAREDE E SUPERFICIE DA CARGA')
disp(' ')
%=========================================================================
Anexo A.1 – Listagem do programa
disp(' ')
disp(' DIMENSOES DA CARGA')
Lc=input(' Informe o comprimento da carga: ');
if (Lc>Li)|(Lc<=0.005)
while (Lc>=Li)|(Lc<=0.005)
if (Lc<=0.005)
errordlg('O valor deve ser maior do que o informado !','DIMENSOES DA CARGA')
disp(' *** O valor deve ser maior do que o informado ! ***')
Lc=input(' Digite novamente o comprimento da carga: ');
else
errordlg(['O comprimento da carga nao pode ser maior que ', num2str(Li),' metro(s) !']...
,'DIMENSOES DA CARGA')
disp([' *** O comprimento da carga nao pode ser maior que ', num2str(Li),...
' metro(s) ! ***'])
Lc=input(' Digite novamente o comprimento da carga: ');
disp(' ')
end
end
end
disp(' ')
Wc=input(' Informe a largura da carga: ');
if (Wc>Wi)|(Wc<=0.005)
while (Wc>=Wi)|(Wc<=0.005)
if (Wc<=0.005)
errordlg('O valor deve ser maior do que o informado !','DIMENSOES DA CARGA')
disp(' *** O valor deve ser maior do que o informado ! ***')
Wc=input(' Digite novamente a largura da carga: ');
else
errordlg(['A largura da carga nao pode ser maior que ', num2str(Wi),' metro(s) !']...
,'DIMENSOES DA CARGA')
disp([' *** A largura da carga nao pode ser maior que ', num2str(Wi)...
,' metro(s) ! ***'])
Wc=input(' Digite novamente a largura da carga: ');
disp(' ')
end
end
end
disp(' ')
Hc=input(' Informe a altura da carga: ');
if (Hc>Hi)|(Hc<=0.005)
while (Hc>=Hi)|(Hc<=0.005)
if (Hc<=0.005)
errordlg('O valor deve ser maior do que o informado !','DIMENSOES DA CARGA')
disp([' *** O valor deve ser maior do que o informado ! ***'])
Hc=input(' Digite novamente a altura da carga: ');
else
errordlg(['A altura da carga nao pode ser maior que ', num2str(Hi),' metro(s) !']...
,'DIMENSOES DA CARGA')
disp([' *** A altura da carga nao pode ser maior que ', num2str(Hi),' metro(s) ! ***'])
Hc=input(' Digite novamente a altura da carga: ');
disp(' ')
end
end
end
Ac12=Wc*Lc;
116
Anexo A.1 – Listagem do programa
117
Ac34=Wc*Hc;
Ac56=Lc*Hc;
Ac=(Ac12+Ac34+Ac56)*2;
disp(' ')
disp([' Area total da carga e ',num2str(Ac)])
Vc=Wc*Hc*Lc;
Ec=Vc/Ac;
EW=Wi-Wc;
EH=Hi-Hc;
EL=Li-Lc;
Epp=((EW+EH+EL)/6);
Epp=Epp/2;
disp(' ')
disp([' O espaçamento medio entre a carga e a parede interna e ', num2str(Epp),' metros'])
disp(' ')
%============================= Condiçoes Internas =============================
qe=input('
Informe a quantidade de calor retirado no evaporador: ');
vsaid=input(' Informe a velocidade do ar na saida do evaporador em m/s: ');
disp(' ')
tarinterno=input(' Informe a temperatura inicial do ar interno: ');
ta=tarinterno;
tarinterno=tarinterno+273.15;
%=========================== CALCULOS DA CARGA ===========================
alfacarga=input(' Informe o valor da difusidade termica da carga: ');
disp(' ')
tempcarga=input(' Informe a temperatura inicial da carga: ');
TK=tempcarga;
kcarga=input(' Informe o valor da condutividade termica da carga: ');
disp(' ')
tempcarga=tempcarga+273.15;
disp(' ')
disp(' ')
if (vsaid==0)
x2=0;
y2=1;
else
x2=1;
y2=0;
end
vinterna=vsaid/2;
tfluc=(tempcarga+tarinterno)/2;
betac=1/tfluc;
%AR INTERNO
Anexo A.1 – Listagem do programa
118
%=========================== Calculo do Ro ar interno ===========================
roint=351.92*tarinterno^(-1.0017);
%======================== Calculo do Cp, K e Alfa ar interno =======================
cpint=(0.103409*10)+(-0.28488708*10^(-3)*tarinterno)+(0.7816818*10^(-6)*tarinterno^2)...
+(-0.4970786*10^(-9)*tarinterno^3)+(0.1077024*10^(-12)*tarinterno^4);
kint=(-2.276501*10^(-3))+(1.2598485*10^(-4)*tarinterno)+(-1.4815235*10^(-7)*tarinterno^2)...
+(1.73550646*10^(-10)*tarinterno^3)+(-1.066657*10^(-13)*tarinterno^4)+(2.47663035*10^(-17)...
*tarinterno^5);
alfaint=(kint/(roint*cpint))/1000;
%CARGA
%============================= Calculo do Ro carga =============================
roc=351.92*tfluc^(-1.0017);
%========================== Calculo do Cp, K e Alfa carga =========================
cpc=(0.103409*10)+(-0.28488708*10^(-3)*tfluc)+(0.7816818*10^(-6)*tfluc^2)+(-0.4970786*10^(-9)...
*tfluc^3)+(0.1077024*10^(-12)*tfluc^4);
kc=(-2.276501*10^(-3))+(1.2598485*10^(-4)*tfluc)+(-1.4815235*10^(-7)*tfluc^2)+(1.73550646*10^(10)...
*tfluc^3)+(-1.066657*10^(-13)*tfluc^4)+(2.47663035*10^(-17)*tfluc^5);
alfac=(kc/(roc*cpc))/1000;
%==================== Calculo da Viscosidade Dinamica e Cinematica ==================
viscosidadedinc=(-9.8601*10^(-1))+(9.080125*10^(-2)*tfluc)+(-1.17635575*10^(-4)*tfluc^2)...
+(1.2349703*10^(-7)*tfluc^3)+(-5.7971299*10^(-11)*tfluc^4);
viscosidadecinc=viscosidadedinc/(roc*1000000);
%============================ Calculo de Prandtl carga ===========================
prandtlc=viscosidadecinc/alfac;
%========================== Calculo de Reynolds carga ===========================
reyc=vinterna*(Wc/2)/viscosidadecinc;
%=========================================================================
if (tempcarga>tarinterno)
tzc=tempcarga-tarinterno;
else
tzc=tarinterno-tempcarga;
end
%=========================== Calculo de Rayleigh carga ==========================
rac=(9.8*betac*tzc*(Hc^3))/(viscosidadecinc*alfac);
%===================== Calculo de Nusselt para a superficie da carga ===================
nuc=((0.037*(reyc^(4/5)))*(prandtlc^(1/3)))*x2...
+(((0.825+((0.387*(rac^(1/6)))/((1+((0.492/prandtlc)^(9/16)))^(8/27))))^2))*y2;
%=========================================================================
disp(' ');
119
Anexo A.1 – Listagem do programa
disp(' Superficie Carga')
disp('
RO
Cp
K
Alfa
disp([roc;cpc;kc;alfac;viscosidadecinc]')
fprintf('\n')
disp('
Prandtl Reynolds
Rayleigh 1
disp([prandtlc;reyc;rac;nuc;]')
fprintf('\n')
V
')
Rayleigh 2
Nusselt ')
%=========================================================================
hc=(nuc*kc)/(y2*Hc+x2*(Wc/2));
disp(['Coeficiente de transferencia de calor por convecçao para superficie da carga :', num2str(hc)]);
disp(' ')
%=========================================================================
%=========================== CALCULOS INTERNOS ===========================
disp(' ')
for c=1:6
disp(' ')
disp('Parede 1: teto ')
disp('Parede 2: piso ')
disp('Parede 3: lateral direita ')
disp('Parede 4: lateral esquerda ')
disp('Parede 5: traseira ')
disp('Parede 6: frente ')
fprintf('\n')
tparedein(c)=input([' Informe a temperatura da superficie interna da parede ', num2str(c)...
,': ']);
tparedein(c)=tparedein(c)+273.15;
%=========================================================================
disp(' ');
if (vinterna==0)
x7=0;
y7=1;
if (c==1)|(c==2)
z7=1;
s7=0;
else
z7=0;
s7=1;
end
else
x7=1;
y7=0;
if (c==1)|(c==2)|(c==5)|(c==6)
z7=1;
s7=0;
else
z7=0;
s7=1;
end
end
Lint(c)=Li*(x7)*(z7)+(Hi)*(x7)*(s7)+((Li*Wi)/((2*Wi)+(2*Li)))*(y7)*(z7)+Hi*(y7)*(s7);
tflui(c)=(tparedein(c)+tarinterno)/2;
betai(c)=1/tflui(c);
Anexo A.1 – Listagem do programa
120
%INTERNO
%============================ Calculo do Ro interno =============================
roi(c)=351.92*tflui(c)^(-1.0017);
%========================= Calculo do Cp, K e Alfa interno =========================
cpi(c)=(0.103409*10)+(-0.28488708*10^(-3)*tflui(c))+(0.7816818*10^(-6)*tflui(c)^2)+(0.4970786*10^(-9)*tflui(c)^3)+(0.1077024*10^(-12)*tflui(c)^4);
ki(c)=(-2.276501*10^(-3))+(1.2598485*10^(-4)*tflui(c))+(-1.4815235*10^(7)*tflui(c)^2)+(1.73550646*10^(-10)*tflui(c)^3)+(-1.066657*10^(-13)*tflui(c)^4)+(2.47663035*10^(17)*tflui(c)^5);
alfai(c)=(ki(c)/(roi(c)*cpi(c)))/1000;
%=================== Calculo da Viscosidade Dinamica e Cinematica ====================
viscosidadedini(c)=(-9.8601*10^(-1))+(9.080125*10^(-2)*tflui(c))+(-1.17635575*10^(-4)*tflui(c)^2)...
+(1.2349703*10^(-7)*tflui(c)^3)+(-5.7971299*10^(-11)*tflui(c)^4);
viscosidadecini(c)=viscosidadedini(c)/(roi(c)*1000000);
%========================== Calculo de Prandtl interno ============================
prandtli(c)=viscosidadecini(c)/alfai(c);
%========================== Calculo de Reynolds interno ==========================
reyi(c)=vinterna*Lint(c)/viscosidadecini(c);
%=========================================================================
if (tparedein(c)>tarinterno)
tz(c)=tparedein(c)-tarinterno;
else
tz(c)=tarinterno-tparedein(c);
end
%============================= Calculo de Rayleigh interno========================
rai(c)=((9.8*betai(c)*tz(c)*(Lint(c)^3))/(viscosidadecini(c)*alfai(c)));
%=========================== Teste da velocidade relativa =========================
if (vinterna==0)
x8=0;
if (c==1)
y8=0;
z8=0;
s8=1;
else if (c==2)
y8=0;
z8=1;
s8=0;
end
if (c==3)|(c==4)|(c==5)|(c==6)
y8=1;
z8=0;
s8=0;
end
end
disp(' ')
Anexo A.1 – Listagem do programa
121
else
x8=1;y8=0;z8=0;s8=0;
end
%============================= Calculo de Nusselt ==============================
nui(c)=((0.037*(reyi(c)^(4/5)))*(prandtli(c)^(1/3)))*x8...
+(((0.825+((0.387*(rai(c)^(1/6)))/((1+((0.492/prandtli(c))^(9/16)))^(8/27))))^2))*y8...
+(0.27*(rai(c)^(1/4)))*s8+(0.15*(rai(c)^(1/3)))*z8;
%=========================================================================
disp(' ');
disp(' Superficie Interna')
disp('
RO
Cp
K
Alfa
V ')
disp([roi(c);cpi(c);ki(c);alfai(c);viscosidadecini(c)]')
fprintf('\n')
disp('
Prandtl Reynolds
Rayleigh
Nusselt ')
disp([prandtli(c);reyi(c);rai(c);nui(c);]')
fprintf('\n')
%=========================================================================
hi(c)=(nui(c)*ki(c))/(Lint(c));
disp(['
Coeficiente interno de convecçao e ', num2str(hi(c))]);
disp(' ')
absortc=
emissivc=
qr(c)=(absortc*tal*(tparedein(c))^4)-(emissivc*tal*(tempcarga)^4);
end
qr= (qr(1)+qr(2)+qr(3)+qr(4)+qr(5)+qr(6))/6;
%============================== EVAPORADOR =============================
if qe>0
kcarga=kint;
alfacarga=alfaint;
end
%=========================================================================
t1=tparedein(1);
t2=tparedein(2);
t3=tparedein(3);
t4=tparedein(4);
t5=tparedein(5);
t6=tparedein(6);
%============================= H da parede interna ==============================
he1=he(1);
he2=he(2);
he3=he(3);
he4=he(4);
he5=he(5);
he6=he(6);
hi1=hi(1);
hi2=hi(2);
hi3=hi(3);
Anexo A.1 – Listagem do programa
122
hi4=hi(4);
hi5=hi(5);
hi6=hi(6);
%==========================Fluxo de calor sobre a carga ===========================
qc=hc*(tarinterno-tempcarga);
qrc=qc+qr;
qri=-(qrc*Ac/6)/Ai12;
%========================= Calculo do Calor Solar Absorvido =======================
sigma=(23.45*sin(2*pi*((284+dia)/365)));
Ws=180*(acos(-tan(latitude*pi/180)*tan(sigma*pi/180)))/pi;
bennetta=input(' Informe o coeficiente de Bennett (a): ');
disp(' ')
bennettb=input(' Informe o coeficiente de Bennett (b): ');
disp(' ')
bennettc=input(' Informe o coeficiente de Bennett (c): ');
disp(' ')
gsc=1367;
insolacao=input(' Informe o numero de insolacao diaria em media mensal: ');
disp(' ')
medinsolac=insolacao/mes;
Nthor=0.1333*(Ws);
disp([' Nthor :', num2str(Nthor)]);
hangular=(15*((hora+.5)-12)); % Hora angular
a=0.409+(0.5016*sin((Ws-60)*pi/180)); % Ws em Graus
disp([' a : ', num2str(a)]);
b=0.6609-(0.4767*sin((Ws-60)*pi/180));
disp([' b : ', num2str(b)]);
Rt=(pi/24)*(a+(b*cos(hangular*pi/180)))*((cos(hangular*pi/180)-cos(Ws*pi/180))/(sin(Ws*pi/180)...
-((pi*Ws/180)*cos(Ws*pi/180))));
disp([' Rt : ', num2str(Rt)]);
Hzero=((24*3600*gsc)/pi)*(1+(0.033*(cos((2*pi*dia)/365))))*(cos(latitude*pi/180)*cos(sigma*pi/180)..
.
*sin(Ws*pi/180)+(Ws*pi/180)*sin(latitude*pi/180)*sin(sigma*pi/180));
disp([' H zero : ', num2str(Hzero)]);
H=Hzero*(bennetta+(bennettb*(medinsolac./Nthor))+(bennettc.*altitude));
disp([' H : ', num2str(H)]);
I=H*Rt;
Radsolmedhor=I;
if (Radsolmedhor<0)
Radsolmedhor=0;
end
disp([' Radiaçao Solar em media horaria ', num2str(I),' j/m^2'])
fprintf('\n')
qs=((Radsolmedhor*Ae12)*absortext1)/3600
%=========================================================================
alfaI=input(' Informe a difusidade termica para o isolamento das paredes: ');
fprintf('\n')
kI=input(' Informe a condutividade termica para o isolamento das paredes: ');
fprintf('\n')
erro=input(' Informe a margem de erro(entre 0 e 1): ');
Anexo A.1 – Listagem do programa
123
if (erro>1)|(erro<=0)
while (erro>1)|(erro<=0)
errordlg('Dado incorreto!','MARGEM DE ERRO')
disp(' Dado incorreto!');
erro=input(' Digite o dado novamente: ');
disp(' ')
end
end
fprintf('\n')
Tm11(1)=Text(1);
Tm12(1)=Text(2);
Tm13(1)=Text(3); %externo
Tm14(1)=Text(4);
Tm15(1)=Text(5);
Tm16(1)=Text(6);
Tm1(1)=tparedein(1);
Tm2(1)=tparedein(2);
Tm3(1)=tparedein(3);
Tm4(1)=tparedein(4); %interno
Tm5(1)=tparedein(5);
Tm6(1)=tparedein(6);
Tmm1(1)=tarinterno;
Tc(1)=tempcarga;
%=========================== Determinando o delta T =============================
c1=((EspI(1)*hi1)/kI);
rI1=1/(2*(c1+1));
dt1=(rI1*EspI(1)^2)/alfaI;
c2=((EspI(2)*hi2)/kI);
rI2=1/(2*(c2+1));
dt2=(rI2*EspI(2)^2)/alfaI;
c3=((EspI(3)*hi3)/kI);
rI3=1/(2*(c3+1));
dt3=(rI3*EspI(3)^2)/alfaI;
c4=((EspI(4)*hi4)/kI);
rI4=1/(2*(c4+1));
dt4=(rI4*EspI(4)^2)/alfaI;
c5=((EspI(5)*hi5)/kI);
rI5=1/(2*(c5+1));
dt5=(rI5*EspI(5)^2)/alfaI;
c6=((EspI(6)*hi6)/kI);
rI6=1/(2*(c6+1));
dt6=(rI6*EspI(6)^2)/alfaI;
rI7=1/(2*(((EspI(1)*hre1)/(kI*Ai12))+((EspI(1)*he1)/kI)+1));
dt7=(rI7*EspI(1)^2)/alfaI;
rI8=1/(2*(((he2*EspI(2))/kI)+1));
dt8=(rI8*EspI(2)^2)/alfaI;
rI9=1/(2*(((EspI(3)*hre3)/(kI*Ai34))+((EspI(3)*he3)/kI)+1));
dt9=(rI9*EspI(3)^2)/alfaI;
rI10=1/(2*(((EspI(4)*hre4)/(kI*Ai34))+((EspI(4)*he4)/kI)+1));
Anexo A.1 – Listagem do programa
124
dt10=(rI10*EspI(4)^2)/alfaI;
rI11=1/(2*(((EspI(5)*hre5)/(kI*Ai56))+((EspI(5)*he5)/kI)+1));
dt11=(rI11*EspI(5)^2)/alfaI;
rI12=1/(2*(((EspI(6)*hre6)/(kI*Ai56))+((EspI(6)*he6)/kI)+1));
dt12=(rI12*EspI(6)^2)/alfaI;
rc=1/((2*Ec*hc)/kcarga);
dt13=(rc*Ec^2)/alfacarga;
ci=Epp/(kint*Ai);
ri=1/(ci*(Ai12*(hi1+hi2)+Ai34*(hi3+hi4)+Ai56*(hi5+hi6)+(hc*Ac)));
dt0=(ri*Epp^2)/alfaint;
disp(' Valor do Delta T:')
v=[dt1 dt2 dt3 dt4 dt5 dt6 dt7 dt8 dt9 dt10 dt11 dt12 dt13 dt0]'
deltaT=min(v)*.7;
fprintf('\n')
disp(deltaT)
ttotal=input(' Informe o tempo total de realizaçao do teste (em minutos): ');
ttotal=ttotal*60;
tempo=deltaT;
t=1;
h=waitbar(0,'Calculando, por favor aguarde...','color',[0.92 0.91 0.82]);
while (tempo<=ttotal)
waitbar(tempo/ttotal,h)
rI1=((alfaI*deltaT)/EspI(1)^2);
rI2=((alfaI*deltaT)/EspI(2)^2);
rI3=((alfaI*deltaT)/EspI(3)^2);
rI4=((alfaI*deltaT)/EspI(4)^2);
rI5=((alfaI*deltaT)/EspI(5)^2);
rI6=((alfaI*deltaT)/EspI(6)^2);
rc=((alfacarga*deltaT)/Ec^2);
ri=((alfaint*deltaT)/Epp^2);
ttt=tempo;
prof=0;
X=prof/(2*((alfacarga*ttt)^0.5));
erfcX=1.5577*(2.718281828^(-0.7182*((X+0.7856)^2))); %Incropera;
%erfX=1-(1.5577*(2.718281828^(-0.7182*((X+0.7856)^2)))); %Bejan;
Y=(prof/(2*((alfacarga*ttt)^0.5)))+((hc*((alfacarga*ttt)^0.5))/kcarga);
erfcY=1.5577*(2.718281828^(-0.7182*((Y+0.7856)^2)));
for i=1:500
%================= Temperatura superficie externa, parede (teto) =======================
Tm11(i+1)=(2*rI1)*Tm1(i)+(1(2*rI1*(((EspI(1)*hre1)/(kI*Ai12))+((EspI(1)*he1)/kI)+1)))*Tm11(i)...
+((2*rI1*EspI(1)*hre1)/(kI*Ai12))*tceu+(2*rI1*(EspI(1)*he1)/kI)*tambiente...
+((2*rI1*EspI(1))/(kI*Ai12))*qs;
teste(7,i)=abs((Tm11(i+1)-Tm11(i))/Tm11(i+1));
if (teste(7,i)<=erro)
inter8=Tm11(i);
in8=i;
valor7=Tm11(i+1);
Anexo A.1 – Listagem do programa
125
va7(t)=Tm11(i+1);
if (t==1)
ext1=Tm11(i+1);
end
end
%==================== Temperatura superficie externa, parede 2 (piso) ==================
Tm12(i+1)=(2*rI2)*Tm2(i)+(1-(2*rI2*(((he2*EspI(2))/kI)+1)))*Tm12(i)+((2*rI2*EspI(2)*he2)/kI)...
*tambiente;
teste(8,i)=abs((Tm12(i+1)-Tm12(i))/Tm12(i+1));
if (teste(8,i)<=erro)
inter9=Tm12(i);
in9=i;
valor8=Tm12(i+1);
va8(t)=Tm12(i+1);
if (t==1)
ext2=Tm12(i+1);
end
end
%=================== Temperatura superficie externa, parede 3(traseira) ==================
Tm13(i+1)=(2*rI3)*Tm3(i)+(1-(2*rI3*(((EspI(3)*hre3)/(kI*Ai34))+((EspI(3)*he3)/kI)+1)))...
*Tm13(i)+ ((2*rI3*EspI(3)*hre3)/(kI*Ai34))*tceu+(2*rI3*(EspI(3)*he3)/kI)*tambiente;
teste(9,i)=abs((Tm13(i+1)-Tm13(i))/Tm13(i+1));
if (teste(9,i)<=erro)
inter10=Tm13(i);
in10=i;
valor9=Tm13(i+1);
va9(t)=Tm13(i+1);
if(t==1)
ext3=Tm13(i+1);
end
end
%=================== Temperatura superficie externa, parede 4 (frontal) ==================
Tm14(i+1)=(2*rI4)*Tm4(i)+(1-(2*rI4*(((EspI(4)*hre4)/(kI*Ai34))+((EspI(4)*he4)/kI)+1)))...
*Tm14(i)+ ((2*rI4*EspI(4)*hre4)/(kI*Ai34))*tceu+((2*rI4*EspI(4)*he4)/kI)*tambiente;
teste(10,i)=abs((Tm14(i+1)-Tm14(i))/Tm14(i+1));
if (teste(10,i)<=erro)
inter11=Tm14(i);
in11=i;
valor10=Tm14(i+1);
va10(t)=Tm14(i+1);
if(t==1)
ext4=Tm14(i+1);
end
end
%================ Temperatura superficie externa, parede 5 (lateral esquerda) ===============
Tm15(i+1)=(2*rI5)*Tm5(i)+(1-(2*rI5*(((EspI(5)*hre5)/(kI*Ai56))+((EspI(5)*he5)/kI)+1)))...
*Tm15(i)+((2*rI5*EspI(5)*hre5)/(kI*Ai56))*tceu+((2*rI5*EspI(5)*he5)/kI)*tambiente;
teste(11,i)=abs((Tm15(i+1)-Tm15(i))/Tm15(i+1));
if (teste(11,i)<=erro)
inter12=Tm15(i);
in12=i;
valor11=Tm15(i+1);
va11(t)=Tm15(i+1);
if (t==1)
Anexo A.1 – Listagem do programa
126
ext5=Tm15(i+1);
end
end
%================ Temperatura superficie externa, parede 6 (lateral direita )===============
Tm16(i+1)=(2*rI6)*Tm6(i)+(1(2*rI6*(((EspI(6)*hre6)/(kI*Ai56))+((EspI(6)*he6)/kI)+1)))*Tm16(i)+((2*rI6*EspI(6)*hre6)/(kI...
*Ai56))*tceu+((2*rI6*EspI(6)*he6)/kI)*tambiente;
teste(12,i)=abs((Tm16(i+1)-Tm16(i))/Tm16(i+1));
if (teste(12,i)<=erro)
inter13=Tm16(i);
in13=i;
valor12=Tm16(i+1);
va12(t)=Tm16(i+1);
if (t==1)
ext6=Tm16(i+1);
end
end
%==================== Temperatura superficie interna, parede 1 (teto) ==================
Tm1(i+1)=(2*rI1*c1)*Tmm1(i)+(1-(2*rI1*(c1+1)))*Tm1(i)+(2*rI1)*Tm11(i+1);
%+2*35*rI1+2*qri*rI1;
teste(1,i)=abs((Tm1(i+1)-Tm1(i))/Tm1(i+1));
if (teste(1,i)<=erro)
inter2=Tm1(i);
in2=i;
valor1=Tm1(i+1);
va1(t)=Tm1(i+1);
if (t==1)
int1=Tm1(i+1);
end
end
%=================== Temperatura superficie interna, parede 2 (piso) ===================
Tm2(i+1)=(2*rI2*c2)*Tmm1(i)+(1-(2*rI2*(c2+1)))*Tm2(i)+(2*rI2)*Tm12(i+1);
%+2*35*rI2+2*qri*rI2;
teste(2,i)=abs((Tm2(i+1)-Tm2(i))/Tm2(i+1));
if (teste(2,i)<=erro)
inter3=Tm2(i);
in3=i;
valor2=Tm2(i+1);
va2(t)=Tm2(i+1);
if (t==1)
int2=Tm2(i+1);
end
end
%================== Temperatura superficie interna, parede 3 (traseira) ===================
Tm3(i+1)=(2*rI3*c3)*Tmm1(i)+(1(2*rI3*(c3+1)))*Tm3(i)+(2*rI3)*Tm13(i+1);
%+2*35*rI3+2*qri*rI3;
teste(3,i)=abs((Tm3(i+1)-Tm3(i))/Tm3(i+1));
if (teste(3,i)<=erro)
inter4=Tm3(i);
in4=i;
valor3=Tm3(i+1);
va3(t)=Tm3(i+1);
if (t==1)
int3=Tm3(i+1);
Anexo A.1 – Listagem do programa
127
end
end
%================== Temperatura superficie interna, parede 4 (frontal) ===================
Tm4(i+1)=(2*rI4*c4)*Tmm1(i)+(1(2*rI4*(c4+1)))*Tm4(i)+(2*rI4)*Tm14(i+1);
%+2*35*rI4+2*qri*rI4;
teste(4,i)=abs((Tm4(i+1)-Tm4(i))/Tm4(i+1));
if (teste(4,i)<=erro)
inter5=Tm4(i);
in5=i;
valor4=Tm4(i+1);
va4(t)=Tm4(i+1);
if (t==1)
int4=Tm4(i+1);
end
end
%================ Temperatura superficie interna, parede 5 (lateral esquerda) ===============
Tm5(i+1)=(2*rI5*c5)*Tmm1(i)+(1(2*rI5*(c5+1)))*Tm5(i)+(2*rI5)*Tm15(i+1);
%+2*35*rI5+2*qri*rI5;
teste(5,i)=abs((Tm5(i+1)-Tm5(i))/Tm5(i+1));
if (teste(5,i)<=erro)
inter6=Tm5(i);
in6=i;
valor5=Tm5(i+1);
va5(t)=Tm5(i+1);
if (t==1)
int5=Tm5(i+1);
end
end
%================ Temperatura superficie interna, parede 6 (lateral direita) =================
Tm6(i+1)=(2*rI6*c6)*Tmm1(i)+(1-(2*rI6*(c6+1)))*Tm6(i)+(2*rI6)*Tm16(i+1);
%+2*35*rI6+2*qri*rI1;
teste(6,i)=abs((Tm6(i+1)-Tm6(i))/Tm6(i+1));
if (teste(6,i)<=erro)
inter7=Tm6(i);
in7=i;
valor6=Tm6(i+1);
va6(t)=Tm6(i+1);
if (t==1)
int6=Tm6(i+1);
end
end
%=========================== Temperatura do ar interno ===========================
Tmm1(i+1)=(1-(ri*ci*(Ai12*(hi1+hi2)+Ai34*(hi3+hi4)+Ai56*(hi5+hi6)+(hc*Ac))))*Tmm1(i)...
+(ri*ci*hi1*Ai12)*Tm1(i+1)+(ri*ci*hi2*Ai12)*Tm2(i+1)+(ri*ci*hi3*Ai34)*Tm3(i+1)...
+(ri*ci*hi4*Ai34)*Tm4(i+1)+(ri*ci*hi5*Ai56)*Tm5(i+1)+(ri*ci*hi6*Ai56)*Tm6(i+1)...
+(ri*ci*hc*Ac)*Tc(i)-(ri*ci*qe);
teste(13,i)=abs((Tmm1(i+1)-Tmm1(i))/Tmm1(i+1));
if (teste(13,i)<=erro)
inter1=Tmm1(i);
in1=i;
valor13=Tmm1(i+1);
va13(t)=Tmm1(i+1);
if (t==1)
ar1=Tmm1(i+1);
Anexo A.1 – Listagem do programa
128
end
end
%======================= Temperatura da superficie da carga =========================
%Tc(i+1)=((erfX+((2.718281828^((hc*prof/kcarga)+((alfacarga*ttt*hc^2)/kcarga^2)))*erfcY))...
*(tempcarga-Tmm1(i+1)))+Tmm1(i+1); % Bejan;
Tc(i+1)=((erfcX-((2.718281828^((hc*prof/kcarga)+((alfacarga*ttt*hc^2)/kcarga^2)))*erfcY))...
*(Tmm1(i+1)-tempcarga))+tempcarga; % Incropera;
%Tc(i+1)=((erfcX-((2.718281828^((hc*prof/kcarga)+((alfacarga*ttt*hc^2)/kcarga^2)))*erfcY))...
*(Tmm1(i+1)-tempcarga))+((2*qrc/kcarga)*((alfacarga*ttt/3.1415)^0.5))+tempcarga;
%Tc(i+1)=((2*qrc/kcarga)*((alfacarga*ttt/3.1415)^0.5))+tempcarga;
%Tc(i+1)=(1-((2*rc*Ec*hc)/kcarga))*Tc(i)+((2*rc*Ec*hc)/kcarga)*Tmm1(i+1);
teste(14,i)=abs((Tc(i+1)-Tc(i))/Tc(i+1));
if (teste(14,i)<=erro)
inter14=Tc(i);
in14=i;
valor14=Tc(i+1);
va14(t)=Tc(i+1);
if (t==1)
carg1=Tc(i+1);
end
end
%=========================================================================
tes=[teste(1,i) teste(2,i) teste(3,i) teste(4,i) teste(5,i) teste(6,i) teste(7,i) teste(8,i)...
teste(9,i) teste(10,i) teste(11,i) teste(12,i) teste(13,i) teste(14,i)]';
fprintf('\n')
%disp(' Max teste')
spi=max(tes);
fprintf('\n')
if (spi<=erro)
z=i;
break
end
i=i+1;
end
%================================TABELAS=================================
%disp([' Numero maximo de iteracoes: ', num2str(z)])
%fprintf('\n')
%i=1:z;
%disp(' i Tm11 Tm12
Tm13
Tm14
Tm15
Tm16
Tc')
%disp([i;teste(7,i);teste(8,i);teste(9,i);teste(10,i);teste(11,i);teste(12,i);teste(14,i)]')
%fprintf('\n')
%disp(' i
Tm1
Tm2
Tm3
Tm4
Tm5
Tm6
Tmm1')
%disp([i;teste(1,i);teste(2,i);teste(3,i);teste(4,i);teste(5,i);teste(6,i);teste(13,i)]')
%================================GRAFICOS=================================
%fprintf('\n')
%disp([' Numero de iterações: ', num2str(in1)])
%fprintf('\n')
%disp(' TEMPERATURA DA CARGA');
%disp([' Temperatura de entrada: ', num2str(tempcarga-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor14-273.15)])
%fprintf('\n')
%disp(' TEMPERATURA DO AR INTERNO');
%disp([' Temperatura de entrada: ', num2str(tarinterno-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor13-273.15)])
Anexo A.1 – Listagem do programa
129
%fprintf('\n')
%disp(' Temperatura externa da parede, face 1');
%disp([' Temperatura de entrada: ', num2str(Text(1)-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor7-273.15)])
%fprintf('\n')
%disp(' Temperatura externa da parede, face 2');
%disp([' Temperatura de entrada: ', num2str(Text(2)-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor8-273.15)])
%fprintf('\n')
%disp(' Temperatura externa da parede, face 3');
%disp([' Temperatura de entrada: ', num2str(Text(3)-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor9-273.15)])
%fprintf('\n')
%disp(' Temperatura externa da parede, face 4');
%disp([' Temperatura de entrada: ', num2str(Text(4)-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor10-273.15)])
%fprintf('\n')
%disp(' Temperatura externa da parede, face 5');
%disp([' Temperatura de entrada: ', num2str(Text(5)-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor11-273.15)])
%fprintf('\n')
%disp(' Temperatura externa da parede, face 6');
%disp([' Temperatura de entrada: ', num2str(Text(6)-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor12-273.15)])
%fprintf('\n')
%disp(' Temperatura do nodo interno, face 1');
%disp([' Temperatura de entrada: ', num2str(tparedein(1)-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor1-273.15)])
%fprintf('\n')
%disp(' Temperatura do nodo interno, face 2');
%disp([' Temperatura de entrada: ', num2str(tparedein(2)-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor2-273.15)])
%fprintf('\n')
%disp(' Temperatura do nodo interno, face 3');
%disp([' Temperatura de entrada: ', num2str(tparedein(3)-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor3-273.15)])
%fprintf('\n')
%disp(' Temperatura do nodo interno, face 4');
%disp([' Temperatura de entrada: ', num2str(tparedein(4)-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor4-273.15)])
%fprintf('\n')
%disp(' Temperatura do nodo interno, face 5');
%disp([' Temperatura de entrada: ', num2str(tparedein(5)-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor5-273.15)])
%fprintf('\n')
%disp(' Temperatura do nodo interno, face 6');
%disp([' Temperatura de entrada: ', num2str(tparedein(6)-273.15)])
%disp([' Temperatura encontrada: ', num2str(valor6-273.15)])
%================================ Acrescimo =================================
Anexo A.1 – Listagem do programa
130
hora2=(tempo/3600)+hora;
if (hora2>=24)
hora2=hora2-24;
end
tambiente2=tmax-(deltat/2)+(deltat/2)*cos((15*((hora2)-14)*pi)/180);
pressaosat2=7*10^(-8)*tambiente2^3+7*10^(-8)*tambiente2^2+tambiente2*6*10^(-5)+0.0005;
pressaov2=pressaosat2*umidader;
pressaovv2=pressaov2*1000;
tdp2=0.0065*pressaovv2^5 - 0.1595*pressaovv2^4 + 1.5601*pressaovv2^3 7.9701*pressaovv2^2+25.739...
*pressaovv2 - 12.192;
tambiente2=tambiente2+273.15;
tceu2=tambiente2*(0.8+(tdp2/250))^0.25;
tempceu2(t)=tceu2-273.15;
%=========================== CALCULOS EXTERNOS ==========================
Text2(1)=valor7;
Text2(2)=valor8;
Text2(3)=valor9;
Text2(4)=valor10;
Text2(5)=valor11;
Text2(6)=valor12;
tfluext2(1)=(Text2(1)+tambiente2)/2;
tfluext2(2)=(Text2(2)+tambiente2)/2;
tfluext2(3)=(Text2(3)+tambiente2)/2;
tfluext2(4)=(Text2(4)+tambiente2)/2;
tfluext2(5)=(Text2(5)+tambiente2)/2;
tfluext2(6)=(Text2(6)+tambiente2)/2;
for e=1:6
betae2(e)=1/tfluext2(e);
ro2(e)=351.92*tfluext2(e)^(-1.0017);
cpexterno2(e)=(0.103409*10)+(-0.28488708*10^(-3)*tfluext2(e))+(0.7816818*10^(6)*tfluext2(e)^2)...
+(-0.4970786*10^(-9)*tfluext2(e)^3)+(0.1077024*10^(-12)*tfluext2(e)^4);
kexterno2(e)=(-2.276501*10^(-3))+(1.2598485*10^(-4)*tfluext2(e))+(-1.4815235*10^(-7)...
*tfluext2(e)^2)+(1.73550646*10^(-10)*tfluext2(e)^3)+(-1.066657*10^(-13)*tfluext2(e)^4)...
+(2.47663035*10^(-17)*tfluext2(e)^5);
alfaexterno2(e)=(kexterno2(e)/(ro2(e)*cpexterno2(e)))/1000;
viscosidadedin2(e)=(-9.8601*10^(-1))+(9.080125*10^(-2)*tfluext2(e))+(-1.17635575*10^(-4)...
*tfluext2(e)^2)+(1.2349703*10^(-7)*tfluext2(e)^3)+(-5.7971299*10^(-11)*tfluext2(e)^4);
viscosidadecin2(e)=viscosidadedin2(e)/(ro2(e)*1000000);
prandtl2(e)=viscosidadecin2(e)/alfaexterno2(e);
reynext2(e)=vrelat*L(e)/viscosidadecin2(e);
if (Text2(e)>tambiente2)
Anexo A.1 – Listagem do programa
131
tx2(e)=Text2(e)-tambiente2;
else
tx2(e)=tambiente2-Text2(e);
end
ra2(e)=((9.8*(betae2(e)*tx2(e))*(L(e)^3))/(viscosidadecin2(e)*alfaexterno2(e)));
if (vrelat==0)
x12=0;
if (e==1)
y12=0;
z12=0;
s12=1;
else if (e==2)
y12=0;
z12=1;
s12=0;
end
if (e==3)|(e==4)|(e==5)|(e==6)
y12=1;
z12=0;
s12=0;
end
end
%=========================== Nusselt Externo 2 ==============================
if (vrelat==0) & (Text2(e)>tambiente2)
x12=0;
y12=0;
if (e==1)
z12=1;
s12=0;
else
z12=0;
s12=1;
end
if (e==2)
z12=0;
s12=1;
else
z12=1;
s12=0;
end
if (e==3)|(e==4)|(e==5)|(e==6)
y12=1;
z12=0;
s12=0;
end
end
else
x12=1;y12=0;z12=0;s12=0;
end
nuext2(e)=((0.037*(reynext2(e)^(4/5)))*(prandtl2(e)^(1/3)))*x12...
+(((0.825+((0.387*(ra2(e)^(1/6)))/((1+((0.492/prandtl2(e))^(9/16)))^(8/27))))^2))*y12...
+(0.27*(ra2(e)^(1/4)))*s12+(0.15*(ra2(e)^(1/3)))*z12;
he7(e)=(nuext2(e)*kexterno2(e))/L(e);
%fprintf('\n')
Anexo A.1 – Listagem do programa
132
%disp([' Superficie externa 2, face : ', num2str(e)])
%disp(' ')
%disp('
RO
Cp
K
Alfa
V ')
%disp([ro2(e);cpexterno2(e);kexterno2(e);alfaexterno2(e);viscosidadecin2(e)]')
%fprintf('\n')
%disp('
Prandtl Reynolds
Rayleigh Nusselt ')
%disp([prandtl2(e);reynext2(e);ra2(e);nuext2(e)]')
%fprintf('\n')
%disp([' Coeficiente externo de convecçao: ',num2str(he2(e))]);
end
he21=he7(1);
he22=he7(2);
he23=he7(3);
he24=he7(4);
he25=he7(5);
he26=he7(6);
hre12=(Ae12)*emissividadeext1*tal*(tceu2^2+Text2(1)^2)*(tceu2+Text2(1));
hre32=(Ae34)*emissividadeext*tal*(tceu2^2+Text2(3)^2)*(tceu2+Text2(3));
hre42=(Ae34)*emissividadeext*tal*(tceu2^2+Text2(4)^2)*(tceu2+Text2(4));
hre52=(Ae56)*emissividadeext*tal*(tceu2^2+Text2(5)^2)*(tceu2+Text2(5));
hre62=(Ae56)*emissividadeext*tal*(tceu2^2+Text2(6)^2)*(tceu2+Text2(6));
%===================================AR INTERNO============================
tarinterno2=valor13;
%========================== Calculo do Ro ar interno =============================
roint2=351.92*tarinterno2^(-1.0017);
%======================= Calculo do Cp, K e Alfa ar interno =========================
cpint2=(0.103409*10)+(-0.2848870*10^(-3)*tarinterno2)+(0.7816818*10^(-6)*tarinterno2^2)...
+(-0.4970786*10^(-9)*tarinterno2^3)+(0.1077024*10^(-12)*tarinterno2^4);
kint2=(-2.276501*10^(-3))+(1.2598485*10^(-4)*tarinterno2)+(-1.4815235*10^(-7)*tarinterno2^2)...
+(1.73550646*10^(-10)*tarinterno2^3)+(-1.066657*10^(-13)*tarinterno2^4)+(2.47663035*10^(-17)...
*tarinterno2^5);
alfaint2=(kint2/(roint2*cpint2))/1000;
%================================== CARGA ================================
Tc2=valor14;
tempcarga2=Tc2-273.15;
kcarga2=kcarga;
alfacarga2=alfacarga;
tfluc2=(Tc2+tarinterno2)/2;
betac2=1/tfluc2;
roc2=351.92*tfluc2.^(-1.0017);
cpc2=(0.103409*10)+(-0.2848870*10^(-3)*tfluc2)+(0.7816818*10^(-6)*tfluc2^2)+(-0.4970786*10^(9)...
*tfluc2^3)+(0.1077024*10^(-12)*tfluc2^4);
Anexo A.1 – Listagem do programa
133
kc2=(-2.276501*10^(-3))+(1.2598485*10^(-4)*tfluc2)+(-1.4815235*10^(-7)*tfluc2^2)...
+(1.73550646*10^(-10)*tfluc2^3)+(-1.066657*10^(-13)*tfluc2^4)+(2.47663035*10^(-17)*tfluc2^5);
alfac2=(kc2/(roc2*cpc2))/1000;
viscosidadedinc2=(-9.8601*10^(-1))+(9.080125*10^(-2)*tfluc2)+(-1.17635575*10^(-4)*tfluc2^2)...
+(1.2349703*10^(-7)*tfluc2^3)+(-5.7971299*10^(-11)*tfluc2^4);
viscosidadecinc2=viscosidadedinc2/(roc2*1000000);
prandtlc2=viscosidadecinc2/alfac2;
reyc2=vinterna*(Wc/2)/viscosidadecinc2;
if (Tc2>tarinterno2)
tzc2=Tc2-tarinterno2;
else
tzc2=tarinterno2-Tc2;
end
rac2=(9.8*betac2*tzc2*(Hc^3))/(viscosidadecinc2*alfac2);
nuc2=((0.037*(reyc2^(4/5)))*(prandtlc2^(1/3)))*x2...
+(((0.825+((0.387*(rac2^(1/6)))/((1+((0.492/prandtlc2)^(9/16)))^(8/27))))^2))*y2;
hc2=(nuc2*kc2)/(y2*Hc+x2*Wc/2);
%============================== Calculos Internos ==============================
tparedein2(1)=valor1;
tparedein2(2)=valor2;
tparedein2(3)=valor3;
tparedein2(4)=valor4;
tparedein2(5)=valor5;
tparedein2(6)=valor6;
tflui2(1)=(tparedein2(1)+tarinterno2)/2;
tflui2(2)=(tparedein2(2)+tarinterno2)/2;
tflui2(3)=(tparedein2(3)+tarinterno2)/2;
tflui2(4)=(tparedein2(4)+tarinterno2)/2;
tflui2(5)=(tparedein2(5)+tarinterno2)/2;
tflui2(6)=(tparedein2(6)+tarinterno2)/2;
for f=1:6
betai2(f)=1/tflui2(f);
roi2(f)=351.92*tflui2(f)^(-1.0017);
cpi2(f)=(0.103409*10)+(-0.28488708*10^(-3)*tflui2(f))+(0.7816818*10^(-6)*tflui2(f)^2)...
+(-0.4970786*10^(-9)*tflui2(f)^3)+(0.1077024*10^(-12)*tflui2(f)^4);
ki2(f)=(-2.276501*10^(-3))+(1.2598485*10^(-4)*tflui2(f))+(-1.4815235*10^(-7)*tflui2(f)^2)...
+(1.73550646*10^(-10)*tflui2(f)^3)+(-1.066657*10^(-13)*tflui2(f)^4)+(2.47663035*10^(-17)...
*tflui2(f)^5);
alfai2(f)=(ki2(f)/(roi2(f)*cpi2(f)))/1000;
viscosidadedini2(f)=(-9.8601*10^(-1))+(9.080125*10^(-2)*tflui2(f))+(-1.17635575*10^(-4)...
*tflui2(f)^2)+(1.2349703*10^(-7)*tflui2(f)^3)+(-5.7971299*10^(-11)*tflui2(f)^4);
viscosidadecini2(f)=viscosidadedini2(f)/(roi2(f)*1000000);
Anexo A.1 – Listagem do programa
prandtli2(f)=viscosidadecini2(f)/alfai2(f);
reyi2(f)=vinterna*Lint(f)/viscosidadecini2(f);
if (tparedein2(f)>tarinterno2)
tz2(f)=tparedein2(f)-tarinterno2;
else
tz2(f)=tarinterno2-tparedein2(f);
end
rai2(f)=(9.8*(betai2(f)*tz2(f))*(Lint(f)^3))/(viscosidadecini2(f)*alfai2(f));
if (vinterna==0)
x82=0;
if (f==1)
y82=0;
z82=0;
s82=1;
else if (f==2)
y82=0;
z82=1;
s82=0;
end
if (f==3)|(f==4)|(f==5)|(f==6)
y82=1;
z82=0;
s82=0;
end
end
disp(' ')
else
x82=1;y82=0;z82=0;s82=0;
end
nui2(f)=((0.037*(reyi2(f)^(4/5)))*(prandtli2(f)^(1/3)))*x82...
+(((0.825+((0.387*(rai2(f)^(1/6)))/((1+((0.492/prandtli2(f))^(9/16)))^(8/27))))^2))*y82...
+(0.27*(rai2(f)^(1/4)))*s82+(0.15*(rai2(f)^(1/3)))*z82;
hi7(f)=(nui2(f)*ki2(f))/(Lint(f));
qr2(f)=(absortc*tal*(tparedein2(f))^4)-(emissivc*tal*(Tc2)^4);
end
qr7= (qr2(1)+qr2(2)+qr2(3)+qr2(4)+qr2(5)+qr2(6))/6;
hi21=hi7(1);
hi22=hi7(2);
hi23=hi7(3);
hi24=hi7(4);
hi25=hi7(5);
hi26=hi7(6);
if qe>0
kcarga2=kint2;
alfacarga2=alfaint2;
end
134
Anexo A.1 – Listagem do programa
135
%==========================Fluxo de calor sobre a carga ===========================
qc7=qc;
qrc2=qc7+qr7;
qri2=-(qrc2*Ac/6)/Ai12;
%qrc=(absortc*tal*(tarinterno2)^4)-(emissivc*tal*(Tc2)^4);
%========================= Calculo do Calor Solar Absorvido =======================
hangular=(15*((hora2+.5)-12)); % Hora angular
Rt=(pi/24)*(a+(b*cos(hangular*pi/180)))*((cos(hangular*pi/180)-cos(Ws*pi/180))/(sin(Ws*pi/180)...
-((pi*Ws/180)*cos(Ws*pi/180))));
I=H*Rt;
Radsolmedhor2=I;
if (Radsolmedhor2<0)
Radsolmedhor2=0;
end
qs2=((Radsolmedhor2*Ae12)*absortext1)/3600;
%============================= Fim do Acressimo =============================
c12=((EspI(1)*hi21)/kI);
rrI12=1/(2*(c12+1));
dt122=(rrI12*EspI(1)^2)/alfaI;
c22=((EspI(2)*hi22)/kI);
rrI22=1/(2*(c22+1));
dt222=(rrI22*EspI(2)^2)/alfaI;
c32=((EspI(3)*hi23)/kI);
rrI32=1/(2*(c32+1));
dt322=(rrI32*EspI(3)^2)/alfaI;
c42=((EspI(4)*hi24)/kI);
rrI42=1/(2*(c42+1));
dt422=(rrI42*EspI(4)^2)/alfaI;
c52=((EspI(5)*hi25)/kI);
rrI52=1/(2*(c52+1));
dt522=(rrI52*EspI(5)^2)/alfaI;
c62=((EspI(6)*hi26)/kI);
rrI62=1/(2*(c62+1));
dt622=(rrI62*EspI(6)^2)/alfaI;
rrI72=1/(2*(((EspI(1)*hre12)/(kI*Ai12))+((EspI(1)*he21)/kI)+1));
dt722=(rrI72*EspI(1)^2)/alfaI;
rrI82=1/(2*(((he22*EspI(2))/kI)+1));
dt822=(rrI82*EspI(2)^2)/alfaI;
rrI92=1/(2*(((EspI(3)*hre32)/(kI*Ai34))+((EspI(3)*he23)/kI)+1));
dt922=(rrI92*EspI(3)^2)/alfaI;
rrI102=1/(2*(((EspI(4)*hre42)/(kI*Ai34))+((EspI(4)*he24)/kI)+1));
dt1022=(rrI102*EspI(4)^2)/alfaI;
Anexo A.1 – Listagem do programa
136
rrI112=1/(2*(((EspI(5)*hre52)/(kI*Ai56))+((EspI(5)*he25)/kI)+1));
dt1122=(rrI112*EspI(5)^2)/alfaI;
rrI122=1/(2*(((EspI(6)*hre62)/(kI*Ai56))+((EspI(6)*he26)/kI)+1));
dt1222=(rrI122*EspI(6)^2)/alfaI;
rrc2=1/((2*Ec*hc2)/kcarga2);
dt1322=(rrc2*Ec^2)/alfacarga2;
ci2=Epp/(kint2*Ai);
rri2=1/(ci2*(Ai12*(hi21+hi22)+Ai34*(hi23+hi24)+Ai56*(hi25+hi26)+(hc2*Ac)));
dt022=(rri2*Epp^2)/alfaint2;
v2=[dt122 dt222 dt322 dt422 dt522 dt622 dt722 dt822 dt922 dt1022 dt1122 dt1222 dt1322 dt022]';
deltaT2=min(v2)*.7;
fprintf('\n')
%============================ Fim do Acrescimo ===============================
%======================== Retornando valores calculados ==========================
tempo=tempo+deltaT2;
deltaT=deltaT2;
Tm11(1)=valor7;
Tm12(1)=valor8;
Tm13(1)=valor9;
Tm14(1)=valor10;
Tm15(1)=valor11;
Tm16(1)=valor12;
Tm1(1)=valor1;
Tm2(1)=valor2;
Tm3(1)=valor3;
Tm4(1)=valor4;
Tm5(1)=valor5;
Tm6(1)=valor6;
Tmm1(1)=valor13;
Tc(1)=valor14;
kcarga=kcarga2;
alfacarga=alfacarga2;
hc=hc2;
kc=kc2;
alfac=alfac2;
qrc=qrc2;
qri=qri2;
ci=ci2;
c1=c12;
c2=c22;
c3=c32;
c4=c42;
c5=c52;
c6=c62;
he1=he21;
he2=he22;
he3=he23;
he4=he24;
he5=he25;
Anexo A.1 – Listagem do programa
137
he6=he26;
hre1=hre12;
hre3=hre32;
hre4=hre42;
hre5=hre52;
hre6=hre62;
qs=qs2;
QS(t)=qs2;
he12(t)=he7(1);
he22(t)=he7(2);
he32(t)=he7(3);
he42(t)=he7(4);
he52(t)=he7(5);
he62(t)=he7(6);
tambiente=tambiente2;
tceu=tceu2;
ambiente(t)=tambiente2;
teceu(t)=tceu2;
rI1e(t)=rI1;
rI2e(t)=rI2;
ra1(t)=ra2(1);
ra2(t)=ra2(2);
ra3(t)=ra2(3);
ra4(t)=ra2(4);
ra5(t)=ra2(5);
ra6(t)=ra2(6);
q1(t)=(kI*Ai12*(valor7-valor1))/EspI(1);
q2(t)=(kI*Ai12*(valor8-valor2))/EspI(2);
q3(t)=(kI*Ai34*(valor9-valor3))/EspI(3);
q4(t)=(kI*Ai34*(valor10-valor4))/EspI(4);
q5(t)=(kI*Ai56*(valor11-valor5))/EspI(5);
q6(t)=(kI*Ai56*(valor12-valor6))/EspI(6);
qtotal(t)=q1(t)+q2(t)+q3(t)+q4(t)+q5(t)+q6(t);
hi1=hi21;
hi2=hi22;
hi3=hi23;
hi4=hi24;
hi5=hi25;
hi6=hi26;
alfaint=alfaint2;
kint=kint2;
qEvap(t)=qEV;
t=t+1;
end
close(h)
fprintf('\n')
para=t-1;
%========================== Fim da devoluçao de valores ==========================
fprintf('\n')
Anexo A.1 – Listagem do programa
disp(' RESULTADOS FINAIS');
fprintf('\n')
disp( 'Radiacao solar em media horaria')
disp(Radsolmedhor2)
fprintf('\n')
disp(' TEMPERATURA DA CARGA');
disp([' Valor da temperatura de entrada: ', num2str(TK)])
disp([' Valor da temperatura no segundo calculo: ', num2str(carg1-273.15)])
disp([' Temperatura encontrada: ', num2str(valor14-273.15)])
fprintf('\n')
disp(' TEMPERATURA DO AR INTERNO');
disp([' Valor inicial de temperatura: ', num2str(ta)])
disp([' Valor da temperatura no segundo calculo: ', num2str(ar1-273.15)])
disp([' Valor final encontrado: ', num2str(valor13-273.15)])
fprintf('\n')
disp(' Temperatura superficie externa, parede 1');
disp([' Valor inicial de temperatura: ', num2str(Text(1)-273.15)])
disp([' Valor da temperatura no segundo calculo: ', num2str(ext1-273.15)])
disp([' Valor final encontrado: ', num2str(valor7-273.15)])
fprintf('\n')
disp(' Temperatura superficie externa, parede 2');
disp([' Valor inicial de temperatura: ', num2str(Text(2)-273.15)])
disp([' Valor da temperatura no segundo calculo: ', num2str(ext2-273.15)])
disp([' Valor final encontrado: ', num2str(valor8-273.15)])
fprintf('\n')
disp(' Temperatura superficie externa, parede 3');
disp([' Valor inicial de temperatura: ', num2str(Text(3)-273.15)])
disp([' Valor da temperatura no segundo calculo: ', num2str(ext3-273.15)])
disp([' Valor final encontrado: ', num2str(valor9-273.15)])
fprintf('\n')
disp(' Temperatura superficie externa, parede 4');
disp([' Valor inicial de temperatura: ', num2str(Text(4)-273.15)])
disp([' Valor da temperatura no segundo calculo: ', num2str(ext4-273.15)])
disp([' Valor final encontrado: ', num2str(valor10-273.15)])
fprintf('\n')
disp(' Temperatura superficie externa, parede 5');
disp([' Valor inicial de temperatura: ', num2str(Text(5)-273.15)])
disp([' Valor da temperatura no segundo calculo: ', num2str(ext5-273.15)])
disp([' Valor final encontrado: ', num2str(valor11-273.15)])
fprintf('\n')
disp(' Temperatura superficie externa, parede 6');
disp([' Valor inicial de temperatura: ', num2str(Text(6)-273.15)])
disp([' Valor da temperatura no segundo calculo: ', num2str(ext6-273.15)])
disp([' Valor final encontrado: ', num2str(valor12-273.15)])
fprintf('\n')
disp(' Temperatura superficie interna, parede 1');
disp([' Valor inicial de temperatura: ', num2str(t1-273.15)])
disp([' Valor da temperatura no segundo calculo: ', num2str(int1-273.15)])
disp([' Valor final encontrado: ', num2str(valor1-273.15)])
fprintf('\n')
disp(' Temperatura superficie interna, parede 2');
138
139
Anexo A.1 – Listagem do programa
disp([' Valor inicial de temperatura: ', num2str(t2-273.15)])
disp([' Valor da temperatura no segundo calculo: ', num2str(int2-273.15)])
disp([' Valor final encontrado: ', num2str(valor2-273.15)])
fprintf('\n')
disp(' Temperatura superficie interna, parede 3');
disp([' Valor inicial de temperatura: ', num2str(t3-273.15)])
disp([' Valor da temperatura no segundo calculo: ', num2str(int3-273.15)])
disp([' Valor final encontrado: ', num2str(valor3-273.15)])
fprintf('\n')
disp(' Temperatura superficie interna, parede 4');
disp([' Valor inicial de temperatura: ', num2str(t4-273.15)])
disp([' Valor da temperatura no segundo calculo: ', num2str(int4-273.15)])
disp([' Valor final encontrado: ', num2str(valor4-273.15)])
fprintf('\n')
disp(' Temperatura superficie interna, parede 5');
disp([' Valor inicial de temperatura: ', num2str(t5-273.15)])
disp([' Valor da temperatura no segundo calculo: ', num2str(int5-273.15)])
disp([' Valor final encontrado: ', num2str(valor5-273.15)])
fprintf('\n')
disp('
disp(['
disp(['
disp(['
Temperatura superficie interna, parede 6');
Valor inicial de temperatura: ', num2str(t6-273.15)])
Valor da temperatura no segundo calculo: ', num2str(int6-273.15)])
Valor final encontrado: ', num2str(valor6-273.15)])
%============================ Impressao de Re-calculo ==========================
%============================ Superficie Interna e Carga =========================
disp(' ');
disp(' Superficie Carga 2')
disp('
RO
Cp
K
Alfa
V
disp([roc2;cpc2;kc2;alfac2;viscosidadecinc2]')
fprintf('\n')
disp('
Prandtl
Reynolds
disp([prandtlc2;reyc2;rac2;nuc2;]')
fprintf('\n')
Rayleigh 1
')
Rayleigh 2
Nusselt ')
disp(['Coeficiente de transferencia de calor por convecçao para superficie da carga :',...
num2str(hc2)]);
fprintf('\n')
disp([' Capacidade calorifica dos gases de exaustao: ', num2str(qEX),' W']);
fprintf('\n')
disp([' Efeito de refrigeraçao disponivel no evaporador: ', num2str(qEV),' W']);
fprintf('\n')
esc=(tempo/t)/60;
t2=t;
tempo=tempo/60;
t=1:para;
clf
h=figure(1);
subplot(2,3,6)
plot(t*esc,va12(t)-273.15)
title ('Tm16 (lat. direita) ')
xlabel('Minutos')
140
Anexo A.1 – Listagem do programa
grid on
axis([0 tempo 20 25])
subplot(2,3,5)
plot(t*esc,va11(t)-273.15)
title ('Tm15 (lat. esquerda)')
xlabel('Minutos')
grid on
axis([0 tempo 20 25])
subplot(2,3,4)
plot(t*esc,va10(t)-273.15)
title ('Tm14 (frontal)')
xlabel('Minutos')
grid on
axis([0 tempo 20 25])
subplot(2,3,3)
plot(t*esc,va9(t)-273.15)
title ('Tm13 (traseira)')
xlabel(' ')
grid on
axis([0 tempo 20 25])
subplot(2,3,2)
plot(t*esc,va8(t)-273.15)
title ('Tm12 (piso)')
xlabel(' ')
grid on
axis([0 tempo 20 25])
subplot(2,3,1)
plot(t*esc,va7(t)-273.15)
title ('Tm11 (teto)')
xlabel(' ')
ylabel('Temperaturas Superficies Externas
grid on
axis([0 tempo 20 50])
h=figure(2);
subplot(2,4,8)
plot(t*esc,va13(t)-273.15)
title ('Tmm1 (ar interno)')
xlabel('Minutos')
grid on
axis([0 tempo 0 20])
subplot(2,4,7)
plot(t*esc,va14(t)-273.15)
title ('Tc (carga)')
xlabel('Minutos')
grid on
axis([0 tempo 0 20])
subplot(2,4,6)
plot(t*esc,va6(t)-273.15)
title ('Tm6 (lat. direita)')
xlabel('Minutos')
grid on
axis([0 tempo 0 20])
subplot(2,4,5)
plot(t*esc,va5(t)-273.15)
title ('Tm5 (lat. esquerda)')
xlabel('Minutos')
grid on
axis([0 tempo 0 20])
')
141
Anexo A.1 – Listagem do programa
subplot(2,4,4)
plot(t*esc,va4(t)-273.15)
title ('Tm4 (frontal)')
xlabel(' ')
grid on
axis([0 tempo 0 20])
subplot(2,4,3)
plot(t*esc,va3(t)-273.15)
title ('Tm3 (traseira)')
xlabel(' ')
grid on
axis([0 tempo 0 20])
subplot(2,4,2)
plot(t*esc,va2(t)-273.15)
title ('Tm2 (piso)')
xlabel(' ')
axis([0 tempo 0 20])
grid on
subplot(2,4,1)
plot(t*esc,va1(t)-273.15)
title ('Tm1 (teto)')
xlabel(' ')
ylabel('Temperaturas Superficies Internas
grid on
axis([0 tempo 0 20])
h=figure(3);
subplot(1,2,1)
plot(t*esc,ambiente(t)-273.15,'r-.',t*esc,teceu(t)-273.15,'b')
legend('Ambiente','Ceu')
title ('Temperatura do Ceu e Ambiente')
xlabel('Minutos')
ylabel('Temperatura')
grid on
axis([0 tempo 15 30])
subplot(1,2,2)
plot(t*esc,QS(t))
axis([0 tempo 0 3500])
title ('Calor Solar')
xlabel('Minutos')
ylabel('W')
grid on
h=figure(4);
plot(t*esc,qtotal(t),'b',t*esc,qEvap(t),'r-.')
axis([0 tempo 0 7000])
legend('Carga Termica','Efeito de Refrigeraçao')
title ('Balanço de Calor')
xlabel('Minutos')
ylabel('Watts')
grid on
')