Projeto_TransCal

•
USP-SC

Gabriel de Oliveira Maia
26/06/2024
Prévia do material em texto
UNIVERSIDADE DE SÃO PAULO
ESCOLA DE ENGENHARIA DE SÃO CARLOS
SEM0550 - TRANSFERÊNCIA DE CALOR E MASSA
Projeto - Análise térmica de uma pista de esqui coberta
localizada em São Carlos
Benito Palma Miele Aniceto - 11232250
Breno Seixas Vieira - 10883828
Catarina Luz M I Barros - 11819636
Daniel Ibiapina Barros Monteiro Martins - 11232500
Gabriel de Oliveira Maia - 11819790
Rafael Rodrigues Garcia - 11819682
Ramon Viana Leôncio - 11819744
São Carlos
Julho de 2023
Sumário
1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1 Escopo do projeto . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Modelo adotado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Determinação dos fatores de forma de condução . . . . . . . . . . . . . . . 5
2.1 Fator de forma da parede virada para o sul . . . . . . . . . . . . . . 5
2.2 Fator de forma entre parede inferior e solo . . . . . . . . . . . . . . 7
3 Simulação de resfriamento da câmara sem neve . . . . . . . . . . . . . . . 9
3.1 Determinação do modelo numérico . . . . . . . . . . . . . . . . . . 10
3.2 Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2.1 Temperaturas . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2.2 Coeficientes de convecção . . . . . . . . . . . . . . . . . . 14
3.2.3 Resistências térmicas em regime permanente . . . . . . . . 15
3.2.4 Potência de refrigeração do evaporador . . . . . . . . . . . 15
3.2.5 Constante de tempo . . . . . . . . . . . . . . . . . . . . . 15
4 Simulação da geração de neve . . . . . . . . . . . . . . . . . . . . . . . . . 16
4.1 Modelo numérico . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4.2 Formação da camada de neve . . . . . . . . . . . . . . . . . . . . . 17
4.3 Temperaturas ao longo do tempo . . . . . . . . . . . . . . . . . . . 17
5 Simulação de controle liga/desliga . . . . . . . . . . . . . . . . . . . . . . . 18
5.1 Gráfico das temperaturas . . . . . . . . . . . . . . . . . . . . . . . . 18
5.2 Intervalo dos ciclos liga/desliga . . . . . . . . . . . . . . . . . . . . 18
5.3 Variação máxima da temperatura do ar no ciclo . . . . . . . . . . . 19
5.4 Variação máxima de temperatura da superf́ıcie lateral do cubo de
gelo durante o liga/desliga [oC] decorrido 8 ciclos. . . . . . . . . . . 19
5.5 Novo intervalo do ciclo liga/desliga em [s] decorridos 8 ciclos liga/-
desliga para constante de tempo igual a 1. . . . . . . . . . . . . . . 19
6 Simulação da ocupação máxima na pista . . . . . . . . . . . . . . . . . . . 20
7 Simulação de falta de energia . . . . . . . . . . . . . . . . . . . . . . . . . 23
8 Simulação de degelo dos evaporadores . . . . . . . . . . . . . . . . . . . . . 25
8.1 Plotar temperaturas . . . . . . . . . . . . . . . . . . . . . . . . . . 25
8.2 Temperatura máxima do ar . . . . . . . . . . . . . . . . . . . . . . 25
9 Implementação Matlab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2
1 Introdução
1.1 Escopo do projeto
O trabalho proposto tem como objetivo fazer a análise e a modelagem térmica de
uma pista de esqui coberta localizada na cidade de São Carlos. Na abordagem proposta, são
considerados diversos elementos que influenciam no comportamento térmico do ambiente,
como, por exemplo, a quantidade de pessoas presentes, a geração de neve, a presença de
teleféricos, além também de elementos estruturais, como paredes e estruturas de aço.
A partir de uma modelagem inicial do problema, é posśıvel realizar várias análises
em relação ao comportamento térmico da pista: simulação de resfriamento da câmara,
simulação da geração de neve, simulação de falta de energia, entre outras situações
relevantes para que o projeto de uma pista de esqui possa ter seu funcionamento ideal em
diversas situações.
1.2 Modelo adotado
Para fazer todas as análises propostas, o primeiro passo é a construção de um
modelo que represente as trocas de calor entre os corpos envolvidos no projeto da pista
de esqui, considerando os três meios de troca de calor: condução, convecção e radiação.
Também, alguns elementos, como os motores, as luzes e as pessoas, serão considerados
como cargas térmicas concentradas.
Assim, o modelo adotado está representado no diagrama abaixo:
3
Figura 1 – Representação das trocas de calor na pista de esqui.
Também, para facilitar a análise, foi feito um modelo que representa o sistema por
meio de resistências térmicas:
Figura 2 – Esquema de resistências térmicas da pista de esqui.
4
2 Determinação dos fatores de forma de condução
Da definição do fator de forma, temos que:
q = Sk∆1−2 (1)
Dos dados do exerćıcio, temos que:
klã = 0, 033W/mK kcim = 0, 35W/mK
Utilizando as espessuras, podemos chegar a keq = 0, 045W/mK. Assim, o fator de
forma pode ser calculado por:
S =
q” · A
k ·∆T
Calculando as áreas médias para cada caso e o fluxo de calor para cada caso por
meio de simulação por elementos finitos, podemos chegar ao resultado para o fator de
forma. Para a simulação, o software escolhido foi o Ansys Engineering Software.
OBS: Foi considerado o encaixe entre uma parede e outra, como demonstrado na
figura:
Figura 3 – Camadas e encaixe das paredes
2.1 Fator de forma da parede virada para o sul
Assumindo que as paredes externas e internas são constantes, podemos simular o
fluxo de calor entre elas para diferentes ∆T .
Primeiramente, foi simulado o fluxo de calor para a condição Tint = −5◦C e
Text = 23◦C, com resultados exibidos nas figuras abaixo:
5
Figura 4 – Fluxo de calor na face interna da parede
Figura 5 – Fluxo de calor na face externa da parede
Em seguida, foi calculado o q”para diferentes ∆T , calculando o fator de forma S
para cada um deles:
q”(W/m2) ∆ T (°C) S (m)
7,26 33 1690,56
8,80 40 1690,53
9,90 45 1690,56
19,81 90 1690,55
41,82 190 1690,56
Tabela 1 – Resultados para o fator de forma da parede sul
Fazendo a média dos resultados obtidos, temos que:
SParSul = 1690, 55 m
6
2.2 Fator de forma entre parede inferior e solo
Para cálculo do fator de forma entre a parede inferior e o solo, foi feito o mesmo
processo, com a exceção de que, em vez de colocar a condição de contorno de temperatura
constante na parede interna, foi colocado um plano inferior representando o solo. Do
exerćıcio, temos que:
ksolo = 1W/mK Tsolo = 20◦C
Figura 6 – Simulação da troca de calor entre parede inferior e chão
Foi feita então a simulação inicialmente com Tint = −5◦C, obtendo os resultados:
Figura 7 – Fluxo de calor total entre chão e parede inferior
7
Figura 8 – Fluxo de calor entre chão e face externa da parede inferior
O q”foi então calculado para diferentes valores de ∆T , obtendo os resultados:
q”(W/m2) Delta T (°C) S (m)
5,42 25 4526,96
15,17 70 4526,96
26,01 120 4526,92
-17,34 -80 4526,93
-39,02 -180 4526,95
Tabela 2 – Resultados para o fator de forma da parede inferior
Fazendo a média dos resultados obtidos, temos que:
SParInf = 4526, 94 m
8
3 Simulação de resfriamento da câmara sem neve
Para a simulação do resfriamento inicial da câmara, será considerado que não
há neve nem a presença de pessoas. Também, o sistema de degelo está desligado, e os
teleféricos estarão com seus motores desligados. Assim, serão analisadas as trocas de calor
representadas na figura 1. Os valores de cada dissipação de calor são:
• qlamp = 10kW , dissipação de calor das lâmpadas, que estão sempre ligadas;
• qvent−evap = 2 · 735.5 · 24 = 35.304kW , calor dos ventiladores dos 24 evaporadores
presentes;
• qmotor−evap = 15 · 735.5 · 24 = 264.78kW , potência total dissipada pelos motores dos
ventiladores do evaporadores;
• qevap = 24 · U ·Alig · (Tevap − Tar,int) = 24 · 10 · (Tevap − Tar,int) = 240(Tevap − Tar,int),
calor trocado entre o ar interno e o evaporador;
• qconv,par,int = Aface(4hpar,int + hteto,int+ hpista)(Tar,int − Tpar,int), calor trocado por
convecção entre a parede interna e o ar interno. Aface é a área de uma única face
interna;
• qcond = Skpar(Tpar,ext−Tpar,int) é o calor trocado por condução entre a parede interna
e a parede externa;
• qconv,par,ext = Aface(4hpar,ext + hteto,ext + hchao)(Tamb − Tpar,ext) é o calor trocado por
convecção entre a parede externa e o ambiente;
• qrad,par,ext = εσ Apar,ext(T
4
amb − T 4
par,ext), calor trocado por radiação entre a parede
externa e o ambiente;
• qrad,par,int = εσ Apar,int(T
4
amb − T 4
par,int), calor trocado por radiação entre a parede
interna e o ambiente;
• qrad,gelo = εσ Agelo(T
4
gelo − T 4
ar,int), calor trocado por radiação entre o cubo de gelo e
o ar interno;
• qconv,gelo = Agelohgelo(Tgelo − Tar,int) é o calor trocado por convecção entre o cubo de
gelo e o ar interno.
9
3.1 Determinação do modelo numérico
Para determinar hcubo,gelo, consideraremos que o cubo de 8 m3, ou seja, de lado 2m,
será composto por uma parede horizontal e 4 paredes verticais. Assim, pode-se utilizar as
correlações emṕıricas para convecção natural em placas verticais e horizontais, e hcubo,gelo
será resultante da associação em pararelo das 5 paredes. Portanto:
• hhor =
Kgelo·0.27Ra
1
4
L
2
;
• hvert =
NuL·Kgelo
2·4 , sendo NuL = 0.68 +
0.670Ra
1/4
L
[1+(0.492/Pr)9/16]
4/9 .
Logo, fazendo a associação: hcubo,gelo = hhor·hver
hhor+hver
.
Para determinar o hpista, utilizamos as correlações de parede plana inclinada. Ou
seja:
hpista =
NuL·Kpista
2·4 , sendo NuL = 0.68+
0.670Ra
1/4
L
[1+(0.492/Pr)9/16]
4/9 e considerando a gravidade
como gcos(15o).
Para determinar os coeficientes de convecção das superf́ıcies interna e externa da
câmara, as paredes serão consideradas também como superf́ıcies planas verticais.
A partir dáı, podemos calcular o coeficiente convectivo para cada superf́ıcie, utili-
zando os resultados desenvolvidos por Churchill e Chu.
Para a modelagem da estrutura da câmara, as paredes foram divididas em dois
elementos de temperaturas uniformes: a parede interna, que continha o revestimento interno
de cimento e metade da massa de lã e a parede externa, que continha o revestimento
externo de cimento e a outra metade da lã.
Tendo sido feitas essas considerações, foi posśıvel aplicar as equações de conservação
de energia para o problema e utilizar métodos de aproximação de derivadas para determinar
a temperatura dos diversos elementos do problema ao longo do tempo. As equações
utilizadas foram:
10
dTar,int
dt
=
qlamp + qvent−evap + qmotor−evap + qevap − qconv,gelo − qconv,par,int
(mcp)ar,int
dTgelo
dt
=
qconv,gelo + qrad,gelo
(mcp)gelo
dTpar,int
dt
=
qconv, par,int − qrad,gelo + qcond
(mcp)par,int
dTpar,ext
dt
=
qrad,par , ext + qconv,par , ext − qcond
(mcp)par,ext
Para cada temperatura, a derivada no tempo pode ser aproximada pela forma:
dT
dt
=
T (t)− T (t−∆t)
∆t
A partir da definição de um valor inicial para as temperaturas, o cálculo das
temperaturas pode ser feito iterativamente:
T (t) = T (t−∆t) + ∆t
dT
dt
A partir dáı, então, é posśıvel fazer a implementação do modelo numérico em
MATLAB.
É importante notar, ainda, que para cada parede, o valor de θ é diferente, além de
também variar com a hora do dia. Assim, segue:
θteto =
π
12
· (hora− 6)
Para 12 ≤ hora ≤ 18 : θoeste =
π
12
· (hora− 12)
Para hora < 12 e hora > 18 : θoeste = 0
Para 6 ≤ hora ≤ 12 : θleste =
π
12
· (hora)
Para hora < 6 e hora > 12 : θleste = 0
θnorte = θsul = 0
Agora, iremos encontrar a expressão para cada troca de calor.
Para o ar interior:
(mcp)ar,int
dTar,int
dt
= qlamp + 24 · qevaporador + qnorte−ar,int + qsul−ar,int + qleste−ar,int + qoeste−ar,int+
+ qteto−ar,int + qpista−ar,int + qaço−ar,int
11
Para a parede interna do teto:
(mcp)teto,int
dTteto,int
dt
= −qteto,int−ar,int + qcond,teto − qteto,int−leste,int − qteto,int−oeste,int−
− qteto,int−sul,int − qteto,int−norte,int − qteto,int−pista − qteto,int−aço − qteto,int−tel
Para a parede externa do teto:
(mcp)teto,ext
dTteto,ext
dt
= −qcond,teto + qrad,ext,teto + qconv,teto,ext
Para a parede interna voltada para Leste:
(mcp)leste,int
dTleste,int
dt
= −qleste,int−ar,int + qcond,leste + qteto,int−leste,int + qoeste,int−leste,int−
− qleste,int−sul,int − qleste,int−norte,int − qleste,int−pista − qleste,int−aço − qleste,int−tel
Para a parede externa voltada para Leste:
(mcp)leste,ext
dTleste,ext
dt
= −qcond,leste + qrad,ext,leste + qconv,leste,ext
Para a parede interna voltada para Oeste:
(mcp)oeste,int
dToeste,int
dt
= −qoeste,int−ar,int + qcond,oeste + qteto,int−oeste,int − qoeste,int−leste,int−
− qoeste,int−sul,int − qoeste,int−norte,int − qoeste,int−pista − qoeste,int−aço − qoeste,int−tel
Para a parede externa voltada para Oeste:
(mcp)oeste,ext
dToeste,ext
dt
= −qcond,oeste + qrad,ext,oeste + qconv,oeste,ext
Para a parede interna voltada para Sul:
(mcp)sul,int
dTsul,int
dt
= −qsul,int−ar,int + qcond,sul + qteto,int−sul,int − qleste,int−sul,int+
+ qoeste,int−sul,int + qnorte,int−sul,int − qsul,int−pista − qsul,int−aço − qsul,int−tel
Para a parede externa voltada para Sul:
(mcp)sul,ext
dTsul,ext
dt
= −qcond,sul + qrad,ext,sul + qconv,sul,ext
12
Para a parede interna voltada para Norte:
(mcp)norte,int
dTnorte,int
dt
= −qnorte,int−ar,int + qcond,norte + qteto,int−norte,int + qleste,int−norte,int+
+ qoeste,int−norte,int − qnorte,int−sul,int − qnorte,int−pista − qnorte,int−aço − qnorte,int−tel
Para a parede externa voltada para Norte:
(mcp)norte,ext
dTnorte,ext
dt
= −qcond,norte + qrad,ext,norte + qconv,norte,ext
Para a parede referente ao chão:
(mcp)pista
dTpista
dt
= −qpista−ar + qteto,int−pista + qleste,int−pista + qoeste,int−pista+
+ qsul,int−pista + qnorte,int−pista − qpista−aço − qpista−tel + qsolo
Para o sensor de temperatura do ar:
dTsensor,ar
dt
=
Tar,int−Ts
τsensor,ar
Para o sensor de temperatura do chão:
dTsensor,chão
dt
=
Tchão−Ts
τsensor,chão
Para o evaporador:
Tevaporador = −20C
Para a estrutura de aço:
(mcp)aço
dTaço
dt
= −qaço + qteto−aço + qleste,aço + qoeste,aço+
+ qsul,aço + qnorte,aço − qaço−pista − qaço−tel
Para o teleférico:
(mcp)tel
dTtel
dt
= −qtel + qteto−tel + qleste,tel + qoeste,tel+
+ qsul,tel + qnorte,tel − qtel−pista + qaço−tel
Para o cubo de gelo e para a neve, consideramos que não há a presença dos dois
componentes inicialmente, uma vez que a pista, de ińıcio, não está refrigerada.
13
3.2 Resultados
3.2.1 Temperaturas
A variação das temperaturas dos elementos da tabela 1 está representada no gráfico
da figura 9. Percebe-se que as temperaturas em regime permanente são: Tar,int = −20oC e
Tpar,int = −20oC.
Figura 9 – Temperaturas em função do tempo
3.2.2 Coeficientes de convecção
Os coeficientes de convecção e suas variações com o tempo estão plotados no gráfico
da figura 10. Os valores de regime permanente são: hpista = 0.8W/(m2K), hpar,int =
0.6W/m2K, hpar,ext = 0.23W/m2K e hcubo,gelo = 0.62W/m2K
14
Figura 10 – Coeficientes de convecção
3.2.3 Resistências térmicas em regime permanente
A resistência térmica pode ser calculada a partir do esquema da figura 2. Assim,
tem-se que a resistência equivalente é de 0.00027oC/W
3.2.4 Potência de refrigeração do evaporador
Após a implementação numérica, foi observado que, em regime permanente, a
potência do evaporador é de 32kW.
3.2.5 Constante de tempo
A constante de tempo τ , é definida como sendo o tempo que se leva para um
instrumento alcançar 63, 2% da resposta de regime permanente. A partir da análise dos
gráficos, chega-se à conclusão de que τ = 0.5e6s = 8333.33min = 138.9h.
15
4 Simulação da geração de neve
Para a geração de neve, será considerado que a transformação da água em neve
ocorrerá imediatamente após a sáıda do canhão. Após essa transformação, a neve então se
deposita no chão, aumentando a altura dacamada de neve gradualmente.
As trocas de calor ocorridas estão esquematizadas na figura 14. A partir dessa
definição inicial do problema, é posśıvel então aprofundar a análise.
Figura 11 – Trocas de calor durante a geração de neve
4.1 Modelo numérico
Para a transformação de água em gelo, o calor é:
qsolid = magua[cp,agua(Tagua − Tgelo) + hls,agua]
E para a transformação de gelo em neve:
qneve = mnevecp,neve(Tgelo − Tneve)
Como toda a água é transformada em gelo, temos que magua = mneve.
A massa de água é dependente da vazão: magua = Q · ρagua ·∆t.
Ou seja, a altura da camada de neve é dada pela expressão:
hneve =
Q · ρagua ·∆t
ρneve,fofaApista
Do balanço de energia:
dTneve
dt
=
qconv,neve + qrad,neve − qcond,neve
(mcp)neve
dTar,int
dt
=
qlamp + qvent−evap + qsolid + qneve + qmotor−evap + qevap − qconv,neve − qconv,par,int − qrad,neve
(mcp)ar,int
16
Novamente, para o processo iterativo, será considerado que:
T (t) = T (t−∆t) + ∆t
dT
dt
4.2 Formação da camada de neve
A partir da equação hneve =
Q·ρagua·∆t
ρneve,fofaApista
, basta substituir o valor desejado de hneve
para então determinar o tempo correspondente.
Portanto, sendo: Apista = 600 · 35 = 21000m2, Q = 3L/s, ρagua = 1kg/L e
ρneve,fofa = 150kg/m3, tem-se que, para formar 10 cm de neve::
∆t =
0.1 · 21000 · 150
1 · 3
= 105000s = 29.17 horas.
4.3 Temperaturas ao longo do tempo
Figura 12 – Variações de temperatura durante a geração de neve
17
5 Simulação de controle liga/desliga
Para esse item foi elaborado um controle liga/desliga, o qual desativa o evaporador
quando um sensor de constante de tempo τ = 5min registra temperatura de −780C
e reativa o dispositivo quando o reservatório retorna a uma temperatura de −720C. A
variação de temperatura do sensor pode ser expressa como:
dTs
dt
=
Tar,int − Ts
τ
Durante o peŕıodo em que o evaporador está desligado, sua temperatura está livre
para variar, seu coeficiente global de transferência de calor diminui, e uma nova equação
deve ser adicionada ao sistema, considerando a massa e calor espećıfico desse elemento:
dTevap
dt
=
−qevap
(mcp)evap
Desse modo, a máquina operou por um total de XX ciclos de liga/desliga exibidos
no diagrama da Figura XX e foram coletados os dados referentes ao último ciclo que serão
exibidos adiante.
5.1 Gráfico das temperaturas
Figura 13 – Variação da temperatura com sensor liga/desliga
5.2 Intervalo dos ciclos liga/desliga
Intervalo do ciclo 8: 848469.10s
18
5.3 Variação máxima da temperatura do ar no ciclo
Temperatura mı́nima do ar interno: 266.14 K
Temperatura máxima do ar interno: 270.15 K
Variação máxima de temperatura do ar no ciclo [oC]: 4.010C
5.4 Variação máxima de temperatura da superf́ıcie lateral do cubo de gelo
durante o liga/desliga [oC] decorrido 8 ciclos.
Temperatura mı́nima do cubo de gelo: 267.53 K
Temperatura máxima do cubo de gelo: 280.64 K
Variação máxima de temperatura da superf́ıcie lateral do cubo de gelo durante o
liga/desliga [oC] decorrido 8 ciclos: 13.11°C
5.5 Novo intervalo do ciclo liga/desliga em [s] decorridos 8 ciclos liga/des-
liga para constante de tempo igual a 1.
Intervalo do ciclo 8 para constante de tempo de 1s: 845712.20 s
19
6 Simulação da ocupação máxima na pista
Considerando que a pista esteja em operação com o controlador de temperatura
ativo como no item 4 e a partir das 9hrs da manhã um evento reúna em média 1000
pessoas sobre a pista para assistir uma competição de esqui, plotamos as temperaturas
ao longo do tempo até às 21hrs. Isso considerando que os canhões de neve serão ligados
durante uma hora a partir das 12hrs.
Figura 14 – Temperatura ao longo do tempo
Para calcular a potência elétrica máxima consumida pelas unidades de refrigeração
durante o peŕıodo em que os canhões de neve estão ligados, podemos seguir os seguintes
passos:
• Calcular a quantidade de calor removida para resfriar a água de 250C para 00C.
• Calcular a quantidade de calor removida durante a solidificação da água a 00C.
• Calcular a quantidade de calor removida para resfriar a água congelada até a
temperatura de equiĺıbrio.
20
• Calcular a quantidade de calor para a ocupação de 1000 pessoas
• Somar as quantidades de calor removidas nos passos anteriores para obter a quanti-
dade total de calor removida.
• Calcular a potência elétrica dividindo a quantidade total de calor removida pelo
tempo de operação dos canhões de neve e da permanência das pessoas
Vamos realizar os cálculos considerando os dados fornecidos:
1) Calor removido para resfriar a água de 25 °C para 00C. Como a vazão de água
fornecida por todos os canhões é de 3 L/s, temos:
Q1 = m · Cp ·∆T = 3 · 4180 · 25 = 313500J/s
2) Calor removido durante a solidificação da água a 0 °C:
Q2 = m · hls = 3 · 333000 = 999000J/s
3) Calor removido para resfriar a água congelada até a temperatura de equiĺıbrio:
Q3 = m · Cp−neve ·∆T = 3 · 2100 · 5 = 31500J/s
4) Quantidade total de calor removida para o canhão:
Qcanhao = Q1 +Q2 +Q3 = 313500 + 999000 + 31500 = 1341000J/s
5) Potência elétrica máxima consumida para os canhões ligados por uma hora (3600
segundos):
Pcanhao =
Qcanhao
t
==
1341000
3600
= 372.5kW
6) Quantidade total de calor removida pela presença das pessoas:
Qpessoas = 100 ·N = 100 · 1000 = 100000J/s
7) Potência elétrica máxima consumida para as pessoas presentes por 12 hora
(43.200 segundos):
21
Ppessoas =
Qpessoas
t
==
100000
43200
= 2.3kW
7) Potência elétrica máxima consumida no total:
Ptotal = 372.5 + 2.3 = 374, 8kW
Portanto, a potência elétrica máxima consumida pelas unidades de refrigeração é
de 374.8 kW.
Para calcular o consumo elétrico em kWh durante o peŕıodo de 12 horas, podemos
utilizar a potência elétrica máxima calculada anteriormente (374.8 kW) e multiplicá-la
pelo tempo de operação correspondente para cada caso.
O consumo elétrico é dado pela fórmula: Consumo elétrico (em kWh) = Potência
(em kW) x Tempo (em horas).
No caso, o tempo de operação dos canhões de neve é de 1 hora. Portanto:
Ccanhoes = 372.5kW · 1h = 372.5kWh
No caso, o tempo de operação para as pessoas é de 12 horas. Portanto:
Cpessoas = 2.3kW · 12h = 27, 78kWh
Assim, o consumo elétrico durante o peŕıodo de 12 horas é de 400.27 kWh.
22
7 Simulação de falta de energia
Para a simulação da falta de energia na pista de esqui, consideraremos as mesmas
condições da simulação liga/desliga (Temperatura do ar interno = 5 °C), porém sem as
trocas de calor causadas pelos aparelhos elétricos, e no intervalo de tempo entre 6:00 e
21:00.
Logo, as trocas de calor consideradas serão: radiação ambiente e convecção natural
nas superf́ıcies externas, condução nas paredes e no teto, radiação e convecção natural
entre as superf́ıcies internas e o ar interno, radiação e convecção natural na pista, condução
na neve, condução no solo, radiação e convecção natural no cubo de gelo e calor gerado
pelas pessoas.
Os objetivos serão determinar a temperatura ao longo do intervalo, bem como os
tempos necessários para a neve e o cubo de gelo começarem a derreter.
O balanço de energia para o ar interno será:
dTar,int
dt
=
qconv,par,int + qconv,pessoas − qconv,pista − qconv,gelo
(mcp)ar,int
Em que
qconv,par,int = Aint,total(4hpar,int + hteto,int)(Tar,int − Tpar,int) =
= 33700(4.0, 6 + 0, 6)(Tar,int − Tpar,int) = 101100(Tar,int − Tpar,int)
qconv,pessoas = 100.100 = 104W
qconv,pista = Apistahpista(Tar,int − Tpista) =
= 21000.0, 8(Tar,int − Tpista) = 16800(Tar,int − Tpista)
qconv,gelo =
5
6
Agelohgelo(Tar,int − Tgelo) =
= 20.0, 62(Tar,int − Tgelo) = 12, 4(Tar,int − Tgelo)
Logo:
dTar,int
dt
=
101100(Tar,int − Tpar,int) + 104 − 16800(Tar,int − Tpista)− 12, 4(Tar,int − Tgelo)
2, 1.108
Obtendo-se a equação:
T (t) =
5
421438
e
t+8400ln(10)
2100
Para a neve, o balanço de energia é:
dTneve
dt
=
qconv,pista+ qcond,solo + qrad,par,int − qcond,gelo + qrad,gelo
(mcp)neve
23
Já para o cubo de gelo, obtemos a seguinte equação:
dTgelo
dt
=
qconv,gelo + qcond,pista + qrad,par,int + qrad,pista
(mcp)gelo
24
8 Simulação de degelo dos evaporadores
Para a simulação do degelo dos evaporadores, consideramos o controle liga/desliga
ativado, e agora desligamos cada evaporador durante uma hora intercaladamente, ligando
a resistência elétrica do evaporador [5kW] por 20 minutos. Assim o UA do evaporador
muda para 200W/ºC (ventiladores desligados). Por fim, também existe um termostato
que desliga a resistência caso a temperatura do evaporador ultrapasse 40ºC.
8.1 Plotar temperaturas
Ao plotar o gráfico das temperaturas no tempo até o sistema completar 5 ciclos,
temos:
Figura 15 – Temperaturas até o sistema completar 5 ciclos
8.2 Temperatura máxima do ar
Durante o degelo, a temperatura máxima do ar encontrada foi de 289,15K, ou
16,36ºC.
25
9 Implementação Matlab
1
2 %QUESTAO 2
3
4 %%%%%% DADOS
5 g = 9.81; %Aceleracao da gravidade
6 sigma = 5.67e-8; %Constante de Stefan -Boltzmann
7
8 S = 1690.55; %Fator de forma da pista
9 E = 0.90; %Emissividade das superficies
10
11 T_amb = 23 + 273.15; %Temperatura do ambiente
12 d_ar = 1.17; %Densidade do ar
13 cp_ar = 1003.805; %Calor especifico do ar
14 k_ar = 0.02; %Condutividade termica do ar
15 visc_din = 16e-6; %Viscosidade dinamica do ar
16 visc_cin = visc_din/d_ar; %Viscosidade cinematica do ar
17 Pr = cp_ar*visc_din/k_ar; %Numero de Prandtl
18 alfa = k_ar/(d_ar*cp_ar); %Difusividade termica do ar
19
20
21 k_la = 0.033; %Condutividade termica da la de rocha
22 d_la = 32; %Densidade da la de rocha
23 cp_la = 840; %Calor especifico da la de rocha
24 k_ci = 0.35; %Densidade da placa cimenticia
25 d_ci = 1700; %Densidade da placa cimenticia
26 cp_ci = 800; %Calor especifico da placa cimenticia
27 V_int = 599.72*34.72*9.72; %Volume interno da camara
28 V_ext = 600*35*10; %Volume externo da camara
29 V_ci_int = 2*599.72*34.72*0.02+2*34.72*9.72*0.02+2*599.72*34.72*0.02; %
Volume do revestimento interno de aco
30 V_ci_ext = 2*600*10*0.02+2*35*10*0.02+2*600*35*0.02; %Volume do
revestimento externo de placa cimenticia
31 m_la = (V_ext - V_int)*d_la; %Massa da la de rocha
32 m_ci_int = V_ci_int*d_ci; %Massa do revestimento interno de placa
cimenticia
33 m_ci_ext = V_ci_ext*d_ci; %Massa do revestimento externo de placa
cimenticia
26
34 C_int = m_ci_int*cp_ci + cp_la*m_la /2; %Capacidade termica da parede
interna
35 C_ext = m_ci_ext*cp_ci + cp_la*m_la /2; %Capacidade termica da parede
externa
36
37 N = 1; %Numero de CUBOS
38 A_rad_cubo = 2*2*5; %area de radiacao dO CUBO DE GELO
39 V_cubo = 8; %Volume do cubo de gelo
40 A_cubo = 2*2*6; %Area dO CUBO DE GELO
41 m_cubo = 7360; %Massa dO CUBO DE GELO
42 cp_cubo = 2100; %Calor especifico dO CUBO DE GELO
43 D_cubo = 1.24; %Diametro da aproximacao dO CUBO para uma esfera
44
45 m_ar = (V_int - V_cubo)*d_ar; %Massa de ar dentro da camara
46 UA = 10000; %Coeficiente global do evaporador ligado
47 T_evap = 273.15 -20; %Temperatura do evaporador ligado
48
49 E_gasta_J = 0; %Energia gasta , em Joules
50
51 ra_cubo = g*( D_cubo ^3)/(alfa*visc_cin); %Numero de Rayleigh para cubo
sem temperatura
52 ra_par_int = g*(9.72^3) /(alfa*visc_cin); %Numero de Rayleigh para parede
interna sem temperatura
53 ra_hor_int = g*(16.41^3) /(alfa*visc_cin); %Numero de Rayleigh para
superficies horizontais internas sem temperaturas
54 ra_par_ext = g*(10^3) /(alfa*visc_cin); %Numero de Rayleigh para parede
externa sem temperatura
55 ra_hor_ext = g*(16.54^3) /(alfa*visc_cin); %Numero de Rayleigh para
superficies horizontais externas sem temperatura
56
57 %%%%% TEMPO
58 dt = 1; %Intervalo de tempo entre iteracoes
59 t_fin = 2e6; %Tempo final de simulacao
60 t = 0:dt:t_fin; %Vetor tempo
61
62
63 %%%%%% VETORES
64 T_ar_int = zeros (1,(1 + t_fin/dt)); %Temperatura do ar interno
65 T_cubo = zeros (1,(1 + t_fin/dt)); %Temperatura do cubo
27
66 T_par_int = zeros (1,(1 + t_fin/dt)); %Temperatura da parede interna
67 T_par_ext = zeros (1,(1 + t_fin/dt)); %Temperatura da parede externa
68 h_cubo = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao ddo cubo
69 h_par_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao da parede
interna
70 h_par_ext = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao da parede
externa
71 h_chao_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do chao
interno
72 s_chao_ext = zeros (1,(1 + t_fin/dt)); %fator de forma do chao externo
73 h_teto_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do teto
interno
74 h_teto_ext = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do teto
externo
75
76
77 %%%%% ITERACOES
78 for i = 1:(1 + t_fin/dt)
79 if i ==1 %Condicoes iniciais
80 T_ar_int(i) = T_amb;
81 T_cubo(i) = 273.15;
82 T_par_int(i) = T_amb;
83 T_par_ext(i) = T_amb;
84 q_evap = 0;
85 q_lamp = 0;
86 else
87 q_lamp = 20; %Calor dissipado pela lampada
88 Tf_cubo = (T_cubo(i-1) + T_ar_int(i-1))/2; %Temperatura de filme
do cubo
89 Tf_par_int = (T_par_int(i-1) + T_ar_int(i-1))/2; %Temperatura de
filme das paredes internas
90 Tf_par_ext = (T_par_ext(i-1) + T_amb)/2; %Temperatura de filme
das paredes externas
91
92 B_cubo = 1/ Tf_cubo; %Beta do cubo
93 B_par_int = 1/ Tf_par_int; %Beta das paredes internas
94 B_par_ext = 1/ Tf_par_ext; %Beta das paredes internas
95
28
96 Ra_cubo = ra_cubo*B_cubo*abs(T_cubo(i-1) - T_ar_int(i-1)); %
Numero de Rayleigh do cubo
97 Ra_par_int = ra_par_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(
i-1)); %Numero de Rayleigh das paredes internas
98 Ra_hor_int = ra_hor_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(
i-1)); %Numero de Rayleigh das superficies horizontais internas
99 Ra_par_ext = ra_par_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb); %
Numero de Rayleigh das paredes externas
100 Ra_hor_ext = ra_hor_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb); %
Numero de Rayleigh das superficies horizontais externas
101
102 Nu_cubo = 2 + 0.589*( Ra_cubo ^(1/4))/((1 + (0.469/ Pr)^(9/16))
^(4/9)); %Numero de Nusselt do cubo
103 Nu_par_int = (0.825 + 0.387*( Ra_par_int ^(1/6))/((1+(0.492/ Pr)
^(9/16))^(8/27)))^2; %Numero de Nusselt das paredes internas
104 Nu_par_ext = (0.825 + 0.387*( Ra_par_ext ^(1/6))/((1+(0.492/ Pr)
^(9/16))^(8/27)))^2; %Numero de Nusselt das paredes externas
105 Nu_teto_int = 0.27*( Ra_hor_int ^(1/4)); %Numero de Nusselt do
teto interno
106 Nu_teto_ext = 0.27*( Ra_hor_ext ^(1/4)); %Numero de Nusselt do
teto externo
107 if Ra_hor_int < 10^7
108 Nu_chao_int = 0.54*( Ra_hor_int ^(1/4)); %Numero de Nusselt do
chao interno para escoamento laminar
109 else
110 Nu_chao_int = 0.15*( Ra_hor_int ^(1/3)); %Numero de Nusselt do
chao interno para escoamento turbulento
111 end
112 if Ra_hor_ext < 10^7
113 Nu_chao_ext = 0.54*( Ra_hor_ext ^(1/4));
114 else
115 Nu_chao_ext = 0.15*( Ra_hor_ext ^(1/4));
116 end
117
118 h_cubo(i) = Nu_cubo*k_ar/D_cubo; %Coeficiente de conveccao do
cubo
119 h_par_int(i) = Nu_par_int*k_ar /9.72; %COeficiente de conveccao
das paredes internas
29
120 h_par_ext(i) = Nu_par_ext*k_ar /10; %Coeficiente de conveccao das
paredes externas
121 h_teto_int(i) = Nu_teto_int*k_ar /16.41; %Coeficiente de
conveccao do teto interno
122 h_teto_ext(i) = Nu_teto_ext*k_ar /16.54; %Coeficiente de
conveccao do teto externo
123 h_chao_int(i) = Nu_chao_int*k_ar /16.41; %Coeficiente de
conveccao do chao interno
124 s_chao_ext(i) = Nu_chao_ext*k_ar /16.54;
125
126 q_evap = UA*( T_evap - T_ar_int(i-1)); %Taxa de calor trocado
pelo ar interno com o evaporador
127 q_conv_cubo = h_cubo(i)*N*A_cubo *( T_ar_int(i-1) - T_cubo(i-1));
%Taxa de calor trocado por conveccao pelo cubo com o ar interno
128 q_rad_cubo = E*sigma*A_rad_cubo *( T_par_int(i-1)^4 - T_cubo(i-1)
^4); %Taxa de calor trocado por radiacao pelo cubo com as paredes
internas
129 q_conv_par_int = (9.72*599.72) *(4* h_par_int(i) + h_teto_int(i)+
h_chao_int(i))*( T_ar_int(i-1) - T_par_int(i-1)); %Taxa de calor
trocado por conveccao pela parede interna com o ar interno
130 q_cond = S*k_la*( T_par_ext(i-1) - T_par_int(i-1)); %Taxa de
calor trocado por conducao pela parede interna com a parede externa
131 q_conv_par_ext = (10*600) *(4* h_par_ext(i) + h_teto_ext(i) +
s_chao_ext(i))*(T_amb - T_par_ext(i-1)); %Taxa de calor trocado por
conveccao pela parede externa com o ambiente externo
132 q_rad_par_ext = E*sigma *(2*600*10+2*35*10+600*35) *(T_amb ^4 -
T_par_ext(i-1)^4); %Taxa de calor trocado por radiacao pela parede
externa com o ambiente externo
133
134 dT_ar_int = (q_lamp + q_evap - q_conv_cubo - q_conv_par_int)/(
m_ar*cp_ar); %Taxa de variacao da temperatura do ar interno
135 dT_cubo = (q_conv_cubo + q_rad_cubo)/(N*m_cubo*cp_cubo); %Taxa
de variacao da temperatura do cubo
136 dT_par_int = (q_conv_par_int - q_rad_cubo + q_cond)/C_int; %Taxa
de variacao da temperatura das paredes internas
137 dT_par_ext = (q_rad_par_ext + q_conv_par_ext - q_cond)/C_ext; %
Taxa de variacao da temperatura das paredes externas
138
30
139 T_ar_int(i) = T_ar_int(i-1) + dt*dT_ar_int; %Temperatura do ar
interno
140 T_cubo(i) = T_cubo(i-1) + dt*dT_cubo; %Temperatura do cubo
141 T_par_int(i) = T_par_int(i-1) + dt*dT_par_int; %Temperatura da
parede interna
142 T_par_ext(i) = T_par_ext(i-1) + dt*dT_par_ext; %Temperatura da
parede externa
143 end
144 if T_ar_int(i) > 194.657 %Corresponde a 95% da variacao de
temperatura para o regime permanente
145 E_gasta_J = E_gasta_J + dt*(2* abs(q_evap) + q_lamp); %Energia
gasta
146 end
147 end
148
149 E_gasta_J
150
151 E_gasta_kWh = E_gasta_J /(3.6e+6) %Energia gasta em kWh
152
153 R_ar_int_amb = (T_amb - T_ar_int(end))/q_cond %Resistencia termica entre
o ar interno e o ambiente
154
155 figure %Plot das temperaturas
156 plot(t, T_ar_int -273.15 , t, T_cubo -273.15 , t, T_par_int -273.15 , t,
T_par_ext -273.15);
157 legend('Ar interno ','Cubo de gelo','Parede interna ','Parede externa ')
158 ylabel('T(C)')
159 xlabel('t(s)')
160 grid on
161
162
163 figure %Plot do coeficiente de conveccao da parede interna
164 plot(t,(4* h_par_int + h_chao_int + h_teto_int)/6,t,(4* h_par_ext +
h_teto_ext)/6,t,h_chao_int ,t,h_cubo)
165 legend('Parede Interna ', 'Parede externa ', 'Chao', 'Cubo de gelo')
166 ylabel('h(W/m^2.K)')
167 xlabel('t(s)')
168 grid on
169
31
170 #################################################################
171
172 %QUESTAO 3
173
174 %%%%%% DADOS
175 g = 9.81; %Aceleracao da gravidade
176 sigma = 5.67e-8; %Constante de Stefan -Boltzmann
177
178 S = 1690.55; %Fator de forma da pista
179 E = 0.90; %Emissividade das superficies
180
181 T_amb = 28 +5*sin(2*pi*6)+ 273.15; %Temperatura do ambiente
182 d_ar = 1.17; %Densidade do ar
183 cp_ar = 1003.805; %Calor especifico do ar
184 k_ar = 0.02; %Condutividade termica do ar
185 visc_din = 16e-6; %Viscosidade dinamica do ar
186 visc_cin = visc_din/d_ar; %Viscosidade cinematica do ar
187 Pr = cp_ar*visc_din/k_ar; %Numero de Prandtl
188 alfa = k_ar/(d_ar*cp_ar); %Difusividade termica do ar
189
190
191 k_la = 0.033; %Condutividade termica da la de rocha
192 d_la = 32; %Densidade da la de rocha
193 cp_la = 840; %Calor especifico da la de rocha
194 k_ci = 0.35; %Densidade da placa cimenticia
195 d_ci = 1700; %Densidade da placa cimenticia
196 cp_ci = 800; %Calor especifico da placa cimenticia
197 V_int = 599.72*34.72*9.72; %Volume interno da camara
198 V_ext = 600*35*10; %Volume externo da camara
199 V_ci_int = 2*599.72*34.72*0.02+2*34.72*9.72*0.02+2*599.72*34.72*0.02; %
Volume do revestimento interno de aco
200 V_ci_ext = 2*600*10*0.02+2*35*10*0.02+2*600*35*0.02; %Volume do
revestimento externo de placa cimenticia
201 m_la = (V_ext - V_int)*d_la; %Massa da la de rocha
202 m_ci_int = V_ci_int*d_ci; %Massa do revestimento interno de placa
cimenticia
203 m_ci_ext = V_ci_ext*d_ci; %Massa do revestimento externo de placa
cimenticia
32
204 C_int = m_ci_int*cp_ci + cp_la*m_la /2; %Capacidade termica da parede
interna
205 C_ext = m_ci_ext*cp_ci + cp_la*m_la /2; %Capacidade termica da parede
externa
206
207 Q = 3; %vazao em L/s
208 d_agua = 1; %1Kg/ L, denisdade da agua
209 k_neve = 0.1;
210 d_neve= 150;
211 cp_agua = 4180;
212 h_solid = 333e3;
213 cp_neve = 2100;
214
215 N = 1; %Numero de CUBOS
216 A_rad_cubo = 2*2*5; %Area de radiacao dO CUBO DE GELO
217 V_cubo = 8; %Volume do cubo de gelo
218 A_cubo = 2*2*6; %Area dO CUBO DE GELO
219 m_cubo = 7360; %Massa dO CUBO DE GELO
220 cp_cubo = 2100; %Calor especifico dO CUBO DE GELO
221 D_cubo = 1.24; %Diametro da aproximacao dO CUBO para uma esfera
222
223 m_ar = (V_int - V_cubo)*d_ar; %Massa de ar dentro da camara
224 UA = 10000; %Coeficiente global do evaporador ligado
225 T_evap = 273.15 -20; %Temperatura do evaporador ligado
226
227
228 ra_cubo = g*( D_cubo ^3)/(alfa*visc_cin); %Numero de Rayleigh para cubo
sem temperatura
229 ra_par_int = g*(9.72^3) /(alfa*visc_cin); %Numero de Rayleigh para parede
interna sem temperatura
230 ra_hor_int = g*(16.41^3) /(alfa*visc_cin); %Numero de Rayleigh para
superficies horizontais internas sem temperaturas
231 ra_par_ext = g*(10^3) /(alfa*visc_cin); %Numero de Rayleigh para parede
externa sem temperatura
232 ra_hor_ext = g*(16.54^3) /(alfa*visc_cin); %Numero de Rayleigh para
superficies horizontais externas sem temperatura
233
234 %%%%% TEMPO
235 dt = 1; %Intervalo de tempo entre iteracoes
33
236 t_fin = 1e6; %Tempo final de simulacao
237 t = 0:dt:t_fin; %Vetor tempo
238
239
240 %%%%%% VETORES
241 T_ar_int = zeros (1,(1 + t_fin/dt)); %Temperatura do ar interno
242 T_cubo = zeros (1,(1 + t_fin/dt)); %Temperatura do cubo
243 T_par_int = zeros (1,(1 + t_fin/dt)); %Temperatura da parede interna
244 T_par_ext = zeros (1,(1 + t_fin/dt)); %Temperatura da parede externa
245 h_cubo = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao ddo cubo
246 h_par_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao da parede
interna
247 h_par_ext = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao da parede
externa
248 h_chao_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do chao
interno
249 s_chao_ext = zeros (1,(1 + t_fin/dt)); %fator de forma do chao externo
250 h_teto_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do teto
interno
251 h_teto_ext = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do teto
externo
252
253
254 %%%%% ITERACOES
255 for i = 1:(1 + t_fin/dt)
256 if i ==1 %Condicoes iniciais
257 T_ar_int(i) = -20+273.15;
258 T_cubo(i) = 273.15;
259 T_par_int(i) = -20+273.15;
260 T_par_ext(i) = -20+273.15;
261 q_evap = -3.0504e+03;
262 q_lamp = 20;
263 m_neve = 0;
264 q_solid = 0;
265 q_neve = 0;
266 else
267 q_lamp = 20; %Calor dissipado pela lampada
268 Tf_cubo = (T_cubo(i-1) + T_ar_int(i-1))/2; %Temperatura de filme
do cubo
34
269 Tf_par_int = (T_par_int(i-1) + T_ar_int(i-1))/2; %Temperatura de
filme das paredes internas
270 Tf_par_ext = (T_par_ext(i-1) + T_amb)/2; %Temperatura de filme
das paredes externa
271
272 B_cubo = 1/ Tf_cubo; %Beta do cubo
273 B_par_int = 1/ Tf_par_int; %Beta das paredes internas
274 B_par_ext = 1/ Tf_par_ext; %Beta das paredes internas
275
276 Ra_cubo = ra_cubo*B_cubo*abs(T_cubo(i-1) - T_ar_int(i-1)); %
Numero de Rayleigh do cubo
277 Ra_par_int = ra_par_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(
i-1)); %Numero de Rayleigh das paredes internas
278 Ra_hor_int = ra_hor_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(
i-1)); %Numero de Rayleigh das superficies horizontais internas
279 Ra_par_ext = ra_par_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb); %
Numero de Rayleigh das paredes externas
280 Ra_hor_ext = ra_hor_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb); %
Numero de Rayleigh das superficies horizontais externas
281
282 Nu_cubo = 2 + 0.589*( Ra_cubo ^(1/4))/((1 + (0.469/ Pr)^(9/16))
^(4/9)); %Numero de Nusselt do cubo
283 Nu_par_int = (0.825 + 0.387*(Ra_par_int ^(1/6))/((1+(0.492/ Pr)
^(9/16))^(8/27)))^2; %Numero de Nusselt das paredes internas
284 Nu_par_ext = (0.825 + 0.387*( Ra_par_ext ^(1/6))/((1+(0.492/ Pr)
^(9/16))^(8/27)))^2; %Numero de Nusselt das paredes externas
285 Nu_teto_int = 0.27*( Ra_hor_int ^(1/4)); %Numero de Nusselt do
teto interno
286 Nu_teto_ext = 0.27*( Ra_hor_ext ^(1/4)); %Numero de Nusselt do
teto externo
287 if Ra_hor_int < 10^7
288 Nu_chao_int = 0.54*( Ra_hor_int ^(1/4)); %Numero de Nusselt do
chao interno para escoamento laminar
289 else
290 Nu_chao_int = 0.15*( Ra_hor_int ^(1/3)); %Numero de Nusselt do
chao interno para escoamento turbulento
291 end
292 if Ra_hor_ext < 10^7
293 Nu_chao_ext = 0.54*( Ra_hor_ext ^(1/4));
35
294 else
295 Nu_chao_ext = 0.15*( Ra_hor_ext ^(1/4));
296 end
297
298 h_cubo(i) = Nu_cubo*k_ar/D_cubo; %Coeficiente de conveccao do
cubo
299 h_par_int(i) = Nu_par_int*k_ar /9.72; %COeficiente de conveccao
das paredes internas
300 h_par_ext(i) = Nu_par_ext*k_ar /10; %Coeficiente de conveccao das
paredes externas
301 h_teto_int(i) = Nu_teto_int*k_ar /16.41; %Coeficiente de
conveccao do teto interno
302 h_teto_ext(i) = Nu_teto_ext*k_ar /16.54; %Coeficiente de
conveccao do teto externo
303 h_chao_int(i) = Nu_chao_int*k_ar /16.41; %Coeficiente de
conveccao do chao interno
304 s_chao_ext(i) = Nu_chao_ext*k_ar /16.54;
305
306 q_evap = UA*( T_evap - T_ar_int(i-1)); %Taxa de calor trocado
pelo ar interno com o evaporador
307 q_conv_cubo = h_cubo(i)*N*A_cubo *( T_ar_int(i-1) - T_cubo(i-1));
%Taxa de calor trocado por conveccao pelo cubo com o ar interno
308 q_rad_cubo = E*sigma*A_rad_cubo *( T_par_int(i-1)^4 - T_cubo(i-1)
^4); %Taxa de calor trocado por radiacao pelo cubo com as paredes
internas
309 q_conv_par_int = (9.72*599.72) *(4* h_par_int(i) + h_teto_int(i) +
h_chao_int(i))*( T_ar_int(i-1) - T_par_int(i-1)); %Taxa de calor
trocado por conveccao pela parede interna com o ar interno
310 q_cond = S*k_la*( T_par_ext(i-1) - T_par_int(i-1)); %Taxa de
calor trocado por conducao pela parede interna com a parede externa
311 q_conv_par_ext = (10*600) *(4* h_par_ext(i) + h_teto_ext(i) +
s_chao_ext(i))*(T_amb - T_par_ext(i-1)); %Taxa de calor trocado por
conveccao pela parede externa com o ambiente externo
312 q_rad_par_ext = E*sigma *(2*600*10+2*35*10+600*35) *(T_amb ^4 -
T_par_ext(i-1)^4); %Taxa de calor trocado por radiacao pela parede
externa com o ambiente externo
313
314 m_neve = m_neve + Q*d_agua*dt;
315 q_solid = Q*d_agua*dt*( cp_agua *25+ h_solid);
36
316 q_neve = Q*d_agua*dt*cp_neve *20;
317
318 dT_ar_int = (-q_solid+q_neve+q_lamp + q_evap - q_conv_cubo -
q_conv_par_int)/(m_ar*cp_ar); %Taxa de variacao da temperatura do ar
interno
319 dT_cubo = (q_conv_cubo + q_rad_cubo)/(N*m_cubo*cp_cubo); %Taxa
de variacao da temperatura do cubo
320 dT_par_int = (q_solid -q_neve+q_conv_par_int - q_rad_cubo +
q_cond)/C_int; %Taxa de variacao da temperatura das paredes internas
321 dT_par_ext = (q_rad_par_ext + q_conv_par_ext - q_cond)/C_ext; %
Taxa de variacao da temperatura das paredes externas
322
323 T_ar_int(i) = T_ar_int(i-1) + dt*dT_ar_int; %Temperatura do ar
interno
324 T_cubo(i) = T_cubo(i-1) + dt*dT_cubo; %Temperatura do cubo
325 T_par_int(i) = T_par_int(i-1) + dt*dT_par_int; %Temperatura da
parede interna
326 T_par_ext(i) = T_par_ext(i-1) + dt*dT_par_ext; %Temperatura da
parede externa
327 end
328 end
329
330
331 figure %Plot das temperaturas
332 plot(t, T_ar_int -273.15 , t, T_cubo -273.15 , t, T_par_int -273.15 , t,
T_par_ext -273.15);
333 legend('Ar interno ','Cubo de gelo','Parede interna ','Parede externa ')
334 ylabel('T(C)')
335 xlabel('t(s)')
336 grid on
337
338
339 #################################################################
340
341 %QUESTAO 4
342
343 %%%%%% DADOS
344 g = 9.81; %Aceleracao da gravidade
345 sigma = 5.67e-8; %Constante de Stefan -Boltzmann
37
346 T_amb = 23 + 273.15; %Temperatura do ambiente
347 d_ar = 1.17; %Densidade do ar
348 cp_ar = 1003.805; %Calor especifico do ar
349 k_ar = 0.02; %Condutividade termica do ar
350 visc_din = 16e-6; %Viscosidade dinamica do ar
351 visc_cin = visc_din/d_ar; %Viscosidade cinematica do ar
352 Pr = cp_ar*visc_din/k_ar; %numero de Prandtl
353 alfa = k_ar/(d_ar*cp_ar); %Difusividade termica do ar
354 S = 2484.76; %Fator de forma da pista
355 E = 0.90; %Emissividade das superficies
356 k_la = 0.033; %Condutividade termica da la de rocha
357 d_la = 32; %Densidade da la de rocha
358 cp_la = 840; %Calor especifico da la de rocha
359 k_ci = 0.35; %Densidade da placa cimenticia
360 d_ci = 1700; %Densidade da placa cimenticia
361 cp_ci = 800; %Calor especifico da placa cimenticia
362 V_int = 599.72*34.72*9.72; %Volume interno da camara
363 V_ext = 600*35*10; %Volume externo da camara
364 V_ci_int = 2*599.72*34.72*0.02+2*34.72*9.72*0.02+2*599.72*34.72*0.02; %
Volume do
365 revestimento interno de aco
366 V_ci_ext = 2*600*10*0.02+2*35*10*0.02+2*600*35*0.02; %Volume do
revestimento
367 externo de placa cimenticia
368 m_la = (V_ext - V_int)*d_la; %Massa da la de rocha
369 m_ci_int = V_ci_int*d_ci; %Massa do revestimento interno de placa
cimenticia
370 m_ci_ext = V_ci_ext*d_ci; %Massa do revestimento externo de placa
cimenticia
371 C_int = m_ci_int*cp_ci + cp_la*m_la /2; %Capacidade termica da parede
interna
372 C_ext = m_ci_ext*cp_ci + cp_la*m_la /2; %Capacidade termica da parede
externa
373 N = 1; %numero de CUBOS
374 A_rad_cubo = 2*2*5; %area de radiacao dO CUBO DE GELO
375 V_cubo = 8; %Volume do cubo de gelo
376 A_cubo = 2*2*6; %area dO CUBO DE GELO
377 m_cubo = 7360; %Massa dO CUBO DE GELO
378 cp_cubo = 2100; %Calor especifico dO CUBO DE GELO
38
379 D_cubo = 1.24; %diametro da aproximacao dO CUBO para uma esfera
380 m_ar = (V_int - V_cubo)*d_ar; %Massa de ar dentro da camara
381 UA = 10000; %Coeficiente global do evaporador ligado
382 % T_evap (1) = 273.15 -20; %Temperatura do evaporador ligado
383 E_gasta_J = 0; %Energia gasta , em Joules
384 ra_cubo = g*( D_cubo ^3)/(alfa*visc_cin); %numero de Rayleigh para cubo
sem temperatura
385 ra_par_int = g*(9.72^3) /(alfa*visc_cin); %numero de Rayleigh para parede
interna sem
386 temperatura
387 ra_hor_int = g*(16.41^3) /(alfa*visc_cin); %numero de Rayleigh para
superficies horizontais
388 internas sem temperaturas
389 ra_par_ext = g*(10^3) /(alfa*visc_cin); %numero de Rayleigh para parede
externa sem
390 temperatura
391 ra_hor_ext = g*(16.54^3) /(alfa*visc_cin); %numero de Rayleigh para
superficies horizontais
392 externas sem temperatura
393 %%%%% TEMPO
394 dt = 0.1; %Intervalo de tempo entre iteracoes
395 t_fin = 300000; %Tempo final de simulacao
396 t = 0:dt:t_fin; %Vetor tempo
397 %%%%% DADOS_3
398 tal = 4*60; %4 min
399 % tal = 1; %1 s
400 T_desl = 273.15 -7; %Temp de desligadmento evaporador ( -7 C )
401 T_liga = 273.15 -3; %Temp de religamento do evaporador ( -3 C )
402 ciclo = 1;
403 T_ar_min = 500;
404 T_ar_max = 0;
405 T_cubo_min = 500;
406 T_cubo_max = 0;
407 UA_desl =200;
408 %%%%%% VETORES
409 % T_ar_int = zeros (1,(1 + t_fin/dt)); %Temperatura do ar interno
410 % T_cubo = zeros (1,(1 + t_fin/dt)); %Temperatura do cubo
411 % T_par_int = zeros (1,(1 + t_fin/dt)); %Temperatura da parede interna
412 % T_par_ext = zeros (1,(1 + t_fin/dt)); %Temperatura da parede externa
39
413 % h_cubo = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao ddo cubo
414 % h_par_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao da
parede interna
415 % h_par_ext = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao da
parede externa
416 % h_chao_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do
chao interno
417 % s_chao_ext = zeros (1,(1 + t_fin/dt)); %fator de forma do chao externo
418 % h_teto_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do
teto interno
419 % h_teto_ext = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do
teto externo
420 %%%%% iteracoes
421 for i = 1:(1 + t_fin/dt)
422 ifi ==1 %condicoes iniciais
423 T_ar_int(i) = -5+273.15;
424 T_cubo(i) = 273.15;
425 T_par_int(i) = T_amb;
426 T_par_ext(i) = T_amb;
427 T_s(i) = T_amb;
428 T_evap(i) = 273.15 -20; %Temperatura do evaporador ligado
429 q_evap = 0;
430 q_lamp = 0;
431 elseif T_s(i-1) >= T_desl %Permanece funcionando ate atingir T de
desligamento
432 q_lamp = 10000; %Calor dissipado pela lampada
433 Tf_cubo = (T_cubo(i-1) + T_ar_int(i-1))/2; %Temperatura de filme do cubo
434 Tf_par_int = (T_par_int(i-1) + T_ar_int(i-1))/2; %Temperatura de filme
das paredes
435 internas
436 Tf_par_ext = (T_par_ext(i-1) + T_amb)/2; %Temperatura de filme das
paredes externas
437 B_cubo = 1/ Tf_cubo; %Beta do cubo
438 B_par_int = 1/ Tf_par_int; %Beta das paredes internas
439 B_par_ext = 1/ Tf_par_ext; %Beta das paredes internas
440 Ra_cubo = ra_cubo*B_cubo*abs(T_cubo(i-1) - T_ar_int(i-1)); %numero de
Rayleigh
441 do cubo
40
442 Ra_par_int = ra_par_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(i-1)); %
numero de
443 Rayleigh das paredes internas
444 Ra_hor_int = ra_hor_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(i-1)); %
numero de
445 Rayleigh das superficies horizontais internas
446 Ra_par_ext = ra_par_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb); %numero
de
447 Rayleigh das paredes externas
448 Ra_hor_ext = ra_hor_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb); %numero
de
449 Rayleigh das superficies horizontais externas
450 Nu_cubo = 2 + 0.589*( Ra_cubo ^(1/4))/((1 + (0.469/ Pr)^(9/16))^(4/9)); %
numero de
451 Nusselt do cubo
452 Nu_par_int = (0.825 + 0.387*( Ra_par_int ^(1/6))/((1+(0.492/ Pr)^(9/16))
^(8/27)))^2;
453 %numero de Nusselt das paredes internas
454 Nu_par_ext = (0.825 + 0.387*( Ra_par_ext ^(1/6))/((1+(0.492/ Pr)^(9/16))
^(8/27)))^2;
455 %numero de Nusselt das paredes externas
456 Nu_teto_int = 0.27*( Ra_hor_int ^(1/4)); %numero de Nusselt do teto
interno
457 Nu_teto_ext = 0.27*( Ra_hor_ext ^(1/4)); %numero de Nusselt do teto
externo
458 if Ra_hor_int < 10^7
459 Nu_chao_int = 0.54*( Ra_hor_int ^(1/4)); %numero de Nusselt do chao
interno para
460 escoamento laminar
461 else
462 Nu_chao_int = 0.15*( Ra_hor_int ^(1/3)); %numero de Nusselt do chao
interno para
463 escoamento turbulento
464 end
465 if Ra_hor_ext < 10^7
466 Nu_chao_ext = 0.54*( Ra_hor_ext ^(1/4));
467 else
468 Nu_chao_ext = 0.15*( Ra_hor_ext ^(1/4));
469 end
41
470 h_cubo(i) = Nu_cubo*k_ar/D_cubo; %Coeficiente de conveccao do cubo
471 h_par_int(i) = Nu_par_int*k_ar /9.72; %COeficiente de conveccao das
paredes internas
472 h_par_ext(i) = Nu_par_ext*k_ar /10; %Coeficiente de conveccao das paredes
externas
473 h_teto_int(i) = Nu_teto_int*k_ar /16.41; %Coeficiente de conveccao do
teto interno
474 h_teto_ext(i) = Nu_teto_ext*k_ar /16.54; %Coeficiente de conveccao do
teto externo
475 h_chao_int(i) = Nu_chao_int*k_ar /16.41; %Coeficiente de conveccao do
chao interno
476 s_chao_ext(i) = Nu_chao_ext*k_ar /16.54;
477 q_evap = UA*( T_evap(i-1) - T_ar_int(i-1)); %Taxa de calor trocado pelo
ar interno com
478 o evaporador
479 q_conv_cubo = h_cubo(i)*N*A_cubo *( T_ar_int(i-1) - T_cubo(i-1)); %Taxa de
calor
480 trocado por conveccao pelo cubo com o ar interno
481 q_rad_cubo = E*sigma*A_rad_cubo *( T_par_int(i-1)^4 - T_cubo(i-1)^4); %
Taxa de
482 calor trocado por radiacao pelo cubo com as paredes internas
483 q_conv_par_int = (9.72*599.72) *(4* h_par_int(i) + h_teto_int(i) +
484 h_chao_int(i))*( T_ar_int(i-1) - T_par_int(i-1)); %Taxa de calor trocado
por conveccao pela
485 parede interna com o ar interno
486 q_cond = S*k_la*( T_par_ext(i-1) - T_par_int(i-1)); %Taxa de calor
trocado por
487 conducao pela parede interna com a parede externa
488 q_conv_par_ext = (10*600) *(4* h_par_ext(i) + h_teto_ext(i) + s_chao_ext(i
))*( T_amb -
489 T_par_ext(i-1)); %Taxa de calor trocado por conveccao pela parede
externa com o ambiente
490 externo
491 q_rad_par_ext = E*sigma *(2*600*10+2*35*10+600*35) *(T_amb ^4 -
492 T_par_ext(i-1)^4); %Taxa de calor trocado por radiacao pela parede
externa com o ambiente
493 externo
494 dT_ar_int = (q_lamp + q_evap - q_conv_cubo - q_conv_par_int)/(m_ar*cp_ar
); %Taxa
42
495 de variacao da temperatura do ar interno
496 dT_cubo = (q_conv_cubo + q_rad_cubo)/(N*m_cubo*cp_cubo); %Taxa de
variacao da
497 temperatura do cubo
498 dT_par_int = (q_conv_par_int - q_rad_cubo + q_cond)/C_int; %Taxa de
variacao da
499 temperatura das paredes internas
500 dT_par_ext = (q_rad_par_ext + q_conv_par_ext - q_cond)/C_ext; %Taxa de
variacao da
501 temperatura das paredes externas
502 T_ar_int(i) = T_ar_int(i-1) + dt*dT_ar_int; %Temperatura do ar interno
503 T_cubo(i) = T_cubo(i-1) + dt*dT_cubo; %Temperatura do cubo
504 T_par_int(i) = T_par_int(i-1) + dt*dT_par_int; %Temperatura da parede
interna
505 T_par_ext(i) = T_par_ext(i-1) + dt*dT_par_ext; %Temperatura da parede
externa
506 T_s(i) = T_s(i-1) + dt*( T_ar_int(i-1) - T_s(i-1))/tal;
507 T_evap(i) = T_evap(i-1);
508 if T_cubo(i)>T_cubo_max
509 T_cubo_max=T_cubo(i);
510 end
511 if T_cubo(i)<T_cubo_min
512 T_cubo_min=T_cubo(i);
513 end
514 else
515 break
516 end
517 end
518 t_desl = i*dt; %Registro do tempo em que desligou na primeira vez
519 while ciclo <= 8 %Contagem de ciclos
520 if ciclo == 8
521 t_inicio = (i+1)*dt;
522 end
523 %Fase de aquecimento ate Ts atingir Tliga = -3C
524 while T_s(i-1) <= T_liga | i<size(t)
525 Tf_cubo = (T_cubo(i-1) + T_ar_int(i-1))/2;
526 Tf_par_int = (T_par_int(i-1) + T_ar_int(i-1))/2;
527 Tf_par_ext = (T_par_ext(i-1) + T_amb)/2;
528 B_cubo = 1/ Tf_cubo;
43
529 B_par_int = 1/ Tf_par_int;
530 B_par_ext = 1/ Tf_par_ext;
531 Ra_cubo = ra_cubo*B_cubo*abs(T_cubo(i-1) - T_ar_int(i-1));
532 Ra_par_int = ra_par_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(i-1));
533 Ra_hor_int = ra_hor_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(i-1));
534 Ra_par_ext = ra_par_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb);
535 Ra_hor_ext = ra_hor_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb);
536 Nu_cubo = 2 + 0.589*( Ra_cubo ^(1/4))/((1 + (0.469/ Pr)^(9/16))^(4/9));
537 Nu_par_int = (0.825 + 0.387*( Ra_par_int ^(1/6))/((1+(0.492/ Pr)^(9/16))
^(8/27)))^2;
538 Nu_par_ext = (0.825 + 0.387*( Ra_par_ext ^(1/6))/((1+(0.492/ Pr)^(9/16))
^(8/27)))^2;
539 Nu_teto_int = 0.27*( Ra_hor_int ^(1/4));
540 Nu_teto_ext = 0.27*( Ra_hor_ext ^(1/4));
541 if Ra_hor_int < 10^7
542 Nu_chao_int = 0.54*( Ra_hor_int ^(1/4));
543 else
544 Nu_chao_int = 0.15*( Ra_hor_int ^(1/3));
545 end
546 if Ra_hor_ext < 10^7
547 Nu_chao_ext = 0.54*( Ra_hor_ext ^(1/4));
548 else
549 Nu_chao_ext = 0.15*( Ra_hor_ext ^(1/4));
550 end
551 h_cubo(i) = Nu_cubo*k_ar/D_cubo;
552 h_par_int(i) = Nu_par_int*k_ar /2.4;
553 h_par_ext(i) = Nu_par_ext*k_ar /2.7;
554 h_teto_int(i) = Nu_teto_int*k_ar /0.6;
555 h_teto_ext(i) = Nu_teto_ext*k_ar /0.675;
556 h_chao_int(i) = Nu_chao_int*k_ar /0.6;
557 s_chao_ext(i) = Nu_chao_ext*k_ar /0.675;
558 %Retiramos o calor do evaporador das equacoes
559 q_evap = UA_desl *( T_evap(i-1) - T_ar_int(i-1)); %Taxa de calor trocado
pelo ar interno
560 com o evaporador
561 q_conv_cubo = h_cubo(i)*N*A_cubo *( T_ar_int(i-1) - T_cubo(i-1)); %Taxa de
calor
562 trocado por conveccao pelo cubo com o ar interno
44
563 q_rad_cubo = E*sigma*A_rad_cubo *( T_par_int(i-1)^4 - T_cubo(i-1)^4); %
Taxa de
564 calor trocado por radiacao pelo cubo com as paredes internas
565 q_conv_par_int = (9.72*599.72) *(4* h_par_int(i) + h_teto_int(i) +
566 h_chao_int(i))*( T_ar_int(i-1) - T_par_int(i-1)); %Taxa de calor trocado
por conveccao pela
567 parede interna com o ar interno
568 q_cond = S*k_la*( T_par_ext(i-1) - T_par_int(i-1)); %Taxa de calor
trocado por
569 conducao pela parede interna com a parede externa
570 q_conv_par_ext = (10*600) *(4* h_par_ext(i) + h_teto_ext(i) + s_chao_ext(i
))*( T_amb -
571 T_par_ext(i-1)); %Taxa de calor trocado por conveccao pela parede
externa com o ambiente
572 externo
573 q_rad_par_ext = E*sigma *(2*600*10+2*35*10+600*35) *(T_amb ^4 -
574 T_par_ext(i-1)^4);%Taxa de calor trocado por radiacao pela parede
externa com o ambiente
575 externo
576 dT_ar_int = (q_lamp - q_conv_cubo - q_conv_par_int)/(m_ar*cp_ar);
577 dT_cubo = (q_conv_cubo + q_rad_cubo)/(N*m_cubo*cp_cubo);
578 dT_par_int = (q_conv_par_int - q_rad_cubo + q_cond)/C_int;
579 dT_par_ext = (q_rad_par_ext + q_conv_par_ext - q_cond)/C_ext;
580 dT_evap = (-1)*q_evap /(70*1000);
581 T_ar_int(i) = T_ar_int(i-1) + dt*dT_ar_int;
582 T_cubo(i) = T_cubo(i-1) + dt*dT_cubo;
583 T_par_int(i) = T_par_int(i-1) + dt*dT_par_int;
584 T_par_ext(i) = T_par_ext(i-1) + dt*dT_par_ext;
585 T_s(i) = T_s(i-1) + dt*( T_ar_int(i-1) - T_s(i-1))/tal;
586 T_evap(i) = T_evap(i-1) + dt*dT_evap;
587 if T_cubo(i)>T_cubo_max
588 T_cubo_max=T_cubo(i);
589 end
590 if T_cubo(i)<T_cubo_min
591 T_cubo_min=T_cubo(i);
592 end
593 if ciclo == 8
594 if T_ar_int(i) >= T_ar_max
595 T_ar_max = T_ar_int(i);
45
596 end
597 if T_ar_int(i) <= T_ar_min
598 T_ar_min = T_ar_int(i);
599 end
600 end
601 i = i+1;
602 end
603 %Fase de resfriamento ate Ts atingir Tdesl = -7C
604 while T_s(i-1) >= T_desl
605 T_evap(i) = T_evap (1);
606 Tf_cubo = (T_cubo(i-1) + T_ar_int(i-1))/2;
607 Tf_par_int = (T_par_int(i-1) + T_ar_int(i-1))/2;
608 Tf_par_ext = (T_par_ext(i-1) + T_amb)/2;
609 B_cubo = 1/ Tf_cubo;
610 B_par_int = 1/ Tf_par_int;
611 B_par_ext = 1/ Tf_par_ext;
612 Ra_cubo = ra_cubo*B_cubo*abs(T_cubo(i-1) - T_ar_int(i-1));
613 Ra_par_int = ra_par_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(i-1));
614 Ra_hor_int = ra_hor_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(i-1));
615 Ra_par_ext = ra_par_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb);
616 Ra_hor_ext = ra_hor_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb);
617 Nu_cubo = 2 + 0.589*( Ra_cubo ^(1/4))/((1 + (0.469/ Pr)^(9/16))^(4/9));
618 Nu_par_int = (0.825 + 0.387*( Ra_par_int ^(1/6))/((1+(0.492/ Pr)^(9/16))
^(8/27)))^2;
619 Nu_par_ext = (0.825 + 0.387*( Ra_par_ext ^(1/6))/((1+(0.492/ Pr)^(9/16))
^(8/27)))^2;
620 Nu_teto_int = 0.27*( Ra_hor_int ^(1/4));
621 Nu_teto_ext = 0.27*( Ra_hor_ext ^(1/4));
622 if Ra_hor_int < 10^7
623 Nu_chao_int = 0.54*( Ra_hor_int ^(1/4));
624 else
625 Nu_chao_int = 0.15*( Ra_hor_int ^(1/3));
626 end
627 if Ra_hor_ext < 10^7
628 Nu_chao_ext = 0.54*( Ra_hor_ext ^(1/4));
629 else
630 Nu_chao_ext = 0.15*( Ra_hor_ext ^(1/4));
631 end
632 h_cubo(i) = Nu_cubo*k_ar/D_cubo;
46
633 h_par_int(i) = Nu_par_int*k_ar /2.4;
634 h_par_ext(i) = Nu_par_ext*k_ar /2.7;
635 h_teto_int(i) = Nu_teto_int*k_ar /0.6;
636 h_teto_ext(i) = Nu_teto_ext*k_ar /0.675;
637 h_chao_int(i) = Nu_chao_int*k_ar /0.6;
638 s_chao_ext(i) = Nu_chao_ext*k_ar /0.675;
639 q_evap = UA*( T_evap(i-1) - T_ar_int(i-1)); %Taxa de calor trocado pelo
ar interno com
640 o evaporador
641 q_conv_cubo = h_cubo(i)*N*A_cubo *( T_ar_int(i-1) - T_cubo(i-1)); %Taxa de
calor
642 trocado por conveccao pelo cubo com o ar interno
643 q_rad_cubo = E*sigma*A_rad_cubo *( T_par_int(i-1)^4 - T_cubo(i-1)^4); %
Taxa de
644 calor trocado por radiacao pelo cubo com as paredes internas
645 q_conv_par_int = (9.72*599.72) *(4* h_par_int(i) + h_teto_int(i) +
646 h_chao_int(i))*( T_ar_int(i-1) - T_par_int(i-1)); %Taxa de calor trocado
por conveccao pela
647 parede interna com o ar interno
648 q_cond = S*k_la*( T_par_ext(i-1) - T_par_int(i-1)); %Taxa de calor
trocado por
649 conducao pela parede interna com a parede externa
650 q_conv_par_ext = (10*600) *(4* h_par_ext(i) + h_teto_ext(i) + s_chao_ext(i
))*( T_amb -
651 T_par_ext(i-1)); %Taxa de calor trocado por conveccao pela parede
externa com o ambiente
652 externo
653 q_rad_par_ext = E*sigma *(2*600*10+2*35*10+600*35) *(T_amb ^4 -
654 T_par_ext(i-1)^4); %Taxa de calor trocado por radiacao pela parede
externa com o ambiente
655 externo
656 dT_ar_int = (q_lamp + q_evap - q_conv_cubo - q_conv_par_int)/(m_ar*cp_ar
);
657 dT_cubo = (q_conv_cubo + q_rad_cubo)/( m_cubo*cp_cubo);
658 dT_par_int = (q_conv_par_int - q_rad_cubo + q_cond)/C_int;
659 dT_par_ext = (q_rad_par_ext + q_conv_par_ext - q_cond)/C_ext;
660 T_ar_int(i) = T_ar_int(i-1) + dt*dT_ar_int;
661 T_cubo(i) = T_cubo(i-1) + dt*dT_cubo;
662 T_par_int(i) = T_par_int(i-1) + dt*dT_par_int;
47
663 T_par_ext(i) = T_par_ext(i-1) + dt*dT_par_ext;
664 T_s(i) = T_s(i-1) + dt*( T_ar_int(i-1) - T_s(i-1))/tal;
665 if T_cubo(i)>T_cubo_max
666 T_cubo_max=T_cubo(i);
667 end
668 if T_cubo(i)<T_cubo_min
669 T_cubo_min=T_cubo(i);
670 end
671 if ciclo == 8
672 if T_ar_int(i) >= T_ar_max
673 T_ar_max = T_ar_int(i);
674 end
675 if T_ar_int(i) <= T_ar_min
676 T_ar_min = T_ar_int(i);
677 end
678 end
679 i = i+1;
680 end
681 if ciclo == 8
682 t_final = i*dt;
683 end
684 ciclo = ciclo + 1; %fim do ciclo 1
685 end
686 intervalo = t_final - t_inicio;
687 fprintf (" Intervalo do ciclo 8: %.2f s\n",intervalo)
688 fprintf (" Temperatura minima do ar interno: %.2f K\n",T_ar_min)
689 fprintf (" Temperatura maxima do ar interno: %.2f K\n",T_ar_max)
690 fprintf (" Temperatura minima do cubo de gelo: %.2f K\n",T_cubo_min)
691 fprintf (" Temperatura maxima do cubo de gelo: %.2f K\n",T_cubo_max)
692 tx = (1:dt:(i-1)*dt+(1-dt)); %redimensionar o vetor tempo e ajustar o
passo
693 figure %Plot das temperaturas
694 plot(tx , T_ar_int , tx , T_cubo , tx , T_par_int , tx , T_par_ext);
695 legend('Ar interno ','Cubo de gelo','Parede interna ','Parede externa ')
696 ylabel('T(K)')
697 xlabel('t(s)')
698 grid on
699
700
48
701 #################################################################
702
703 %QUESTAO 5
704
705 %%%%%% DADOS
706 g = 9.81; %Aceleracao da gravidade
707 sigma = 5.67e-8; %Constante de Stefan -Boltzmann
708
709 S = 1690.55; %Fator de forma da pista
710 E = 0.90; %Emissividade das superficies
711
712 T_amb = 28 +5*sin(2*pi*6)+ 273.15; %Temperatura do ambiente
713 d_ar = 1.17; %Densidade do ar
714 cp_ar = 1003.805; %Calor especifico do ar
715 k_ar = 0.02; %Condutividade termica do ar
716 visc_din = 16e-6; %Viscosidade dinamica do ar
717 visc_cin = visc_din/d_ar; %Viscosidade cinematica do ar
718 Pr = cp_ar*visc_din/k_ar; %Numero de Prandtl
719 alfa = k_ar/(d_ar*cp_ar); %Difusividade termica do ar
720
721
722 k_la = 0.033; %Condutividade termica da la de rocha
723 d_la = 32; %Densidade da la de rocha
724 cp_la = 840; %Calor especifico da la de rocha
725 k_ci = 0.35; %Densidade da placa cimenticia
726 d_ci = 1700; %Densidade da placa cimenticia
727 cp_ci = 800; %Calor especifico da placa cimenticia
728 V_int = 599.72*34.72*9.72; %Volume interno da camara
729 V_ext = 600*35*10; %Volume externo da camara
730 V_ci_int = 2*599.72*34.72*0.02+2*34.72*9.72*0.02+2*599.72*34.72*0.02; %
Volume do revestimento interno de aco
731 V_ci_ext = 2*600*10*0.02+2*35*10*0.02+2*600*35*0.02; %Volume do
revestimento externo de placa cimenticia
732 m_la = (V_ext - V_int)*d_la; %Massa da la de rocha
733 m_ci_int = V_ci_int*d_ci; %Massa do revestimento interno de placa
cimenticia
734 m_ci_ext = V_ci_ext*d_ci; %Massa do revestimento externo de placa
cimenticia
49
735 C_int = m_ci_int*cp_ci + cp_la*m_la /2; %Capacidade termica da parede
interna
736 C_ext = m_ci_ext*cp_ci + cp_la*m_la /2; %Capacidade termica da parede
externa
737
738 Q = 3; %vazao em L/s
739 d_agua = 1; %1Kg/ L, denisdade da agua
740 k_neve = 0.1;
741 d_neve= 150;
742 cp_agua = 4180;
743 h_solid = 333e3;
744 cp_neve = 2100;
745
746 N = 1; %Numero de CUBOS
747 A_rad_cubo = 2*2*5; %Area de radiacao dO CUBO DE GELO
748 V_cubo = 8; %Volume do cubo de gelo
749 A_cubo = 2*2*6; %Area dO CUBO DE GELO
750 m_cubo = 7360; %Massa dO CUBO DE GELO
751 cp_cubo = 2100; %Calor especifico dO CUBO DE GELO
752 D_cubo = 1.24; %Diametro da aproximacao dO CUBO para uma esfera
753
754 m_ar = (V_int - V_cubo)*d_ar; %Massa de ar dentro da camara
755 UA = 10000; %Coeficiente global do evaporador ligado
756 T_evap = 273.15 -20; %Temperatura do evaporador ligado
757
758
759 ra_cubo = g*( D_cubo ^3)/(alfa*visc_cin); %Numero de Rayleigh para cubo
sem temperatura
760 ra_par_int= g*(9.72^3) /(alfa*visc_cin); %Numero de Rayleigh para parede
interna sem temperatura
761 ra_hor_int = g*(16.41^3) /(alfa*visc_cin); %Numero de Rayleigh para
superficies horizontais internas sem temperaturas
762 ra_par_ext = g*(10^3) /(alfa*visc_cin); %Numero de Rayleigh para parede
externa sem temperatura
763 ra_hor_ext = g*(16.54^3) /(alfa*visc_cin); %Numero de Rayleigh para
superficies horizontais externas sem temperatura
764
765 %%%%% TEMPO
766 dt = 1; %Intervalo de tempo entre iteracoes
50
767 t_fin = 1e6; %Tempo final de simulacao
768 t = 0:dt:t_fin; %Vetor tempo
769
770
771 %%%%%% VETORES
772 T_ar_int = zeros (1,(1 + t_fin/dt)); %Temperatura do ar interno
773 T_cubo = zeros (1,(1 + t_fin/dt)); %Temperatura do cubo
774 T_par_int = zeros (1,(1 + t_fin/dt)); %Temperatura da parede interna
775 T_par_ext = zeros (1,(1 + t_fin/dt)); %Temperatura da parede externa
776 h_cubo = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao ddo cubo
777 h_par_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao da parede
interna
778 h_par_ext = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao da parede
externa
779 h_chao_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do chao
interno
780 s_chao_ext = zeros (1,(1 + t_fin/dt)); %fator de forma do chao externo
781 h_teto_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do teto
interno
782 h_teto_ext = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do teto
externo
783
784
785 %%%%% ITERACOES
786 for i = 1:(1 + t_fin/dt)
787 if i ==1 %Condicoes iniciais
788 T_ar_int(i) = -20+273.15;
789 T_cubo(i) = 273.15;
790 T_par_int(i) = -20+273.15;
791 T_par_ext(i) = -20+273.15;
792 q_evap = -3.0504e+03;
793 q_lamp = 20;
794 m_neve = 0;
795 q_solid = 0;
796 q_neve = 0;
797 else
798 q_lamp = 20; %Calor dissipado pela lampada
799 Tf_cubo = (T_cubo(i-1) + T_ar_int(i-1))/2; %Temperatura de filme
do cubo
51
800 Tf_par_int = (T_par_int(i-1) + T_ar_int(i-1))/2; %Temperatura de
filme das paredes internas
801 Tf_par_ext = (T_par_ext(i-1) + T_amb)/2; %Temperatura de filme
das paredes externa
802
803 B_cubo = 1/ Tf_cubo; %Beta do cubo
804 B_par_int = 1/ Tf_par_int; %Beta das paredes internas
805 B_par_ext = 1/ Tf_par_ext; %Beta das paredes internas
806
807 Ra_cubo = ra_cubo*B_cubo*abs(T_cubo(i-1) - T_ar_int(i-1)); %
Numero de Rayleigh do cubo
808 Ra_par_int = ra_par_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(
i-1)); %Numero de Rayleigh das paredes internas
809 Ra_hor_int = ra_hor_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(
i-1)); %Numero de Rayleigh das superficies horizontais internas
810 Ra_par_ext = ra_par_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb); %
Numero de Rayleigh das paredes externas
811 Ra_hor_ext = ra_hor_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb); %
Numero de Rayleigh das superficies horizontais externas
812
813 Nu_cubo = 2 + 0.589*( Ra_cubo ^(1/4))/((1 + (0.469/ Pr)^(9/16))
^(4/9)); %Numero de Nusselt do cubo
814 Nu_par_int = (0.825 + 0.387*( Ra_par_int ^(1/6))/((1+(0.492/ Pr)
^(9/16))^(8/27)))^2; %Numero de Nusselt das paredes internas
815 Nu_par_ext = (0.825 + 0.387*( Ra_par_ext ^(1/6))/((1+(0.492/ Pr)
^(9/16))^(8/27)))^2; %Numero de Nusselt das paredes externas
816 Nu_teto_int = 0.27*( Ra_hor_int ^(1/4)); %Numero de Nusselt do
teto interno
817 Nu_teto_ext = 0.27*( Ra_hor_ext ^(1/4)); %Numero de Nusselt do
teto externo
818 if Ra_hor_int < 10^7
819 Nu_chao_int = 0.54*( Ra_hor_int ^(1/4)); %Numero de Nusselt do
chao interno para escoamento laminar
820 else
821 Nu_chao_int = 0.15*( Ra_hor_int ^(1/3)); %Numero de Nusselt do
chao interno para escoamento turbulento
822 end
823 if Ra_hor_ext < 10^7
824 Nu_chao_ext = 0.54*( Ra_hor_ext ^(1/4));
52
825 else
826 Nu_chao_ext = 0.15*( Ra_hor_ext ^(1/4));
827 end
828
829 h_cubo(i) = Nu_cubo*k_ar/D_cubo; %Coeficiente de conveccao do
cubo
830 h_par_int(i) = Nu_par_int*k_ar /9.72; %COeficiente de conveccao
das paredes internas
831 h_par_ext(i) = Nu_par_ext*k_ar /10; %Coeficiente de conveccao das
paredes externas
832 h_teto_int(i) = Nu_teto_int*k_ar /16.41; %Coeficiente de
conveccao do teto interno
833 h_teto_ext(i) = Nu_teto_ext*k_ar /16.54; %Coeficiente de
conveccao do teto externo
834 h_chao_int(i) = Nu_chao_int*k_ar /16.41; %Coeficiente de
conveccao do chao interno
835 s_chao_ext(i) = Nu_chao_ext*k_ar /16.54;
836
837 q_evap = UA*( T_evap - T_ar_int(i-1)); %Taxa de calor trocado
pelo ar interno com o evaporador
838 q_conv_cubo = h_cubo(i)*N*A_cubo *( T_ar_int(i-1) - T_cubo(i-1));
%Taxa de calor trocado por conveccao pelo cubo com o ar interno
839 q_rad_cubo = E*sigma*A_rad_cubo *( T_par_int(i-1)^4 - T_cubo(i-1)
^4); %Taxa de calor trocado por radiacao pelo cubo com as paredes
internas
840 q_conv_par_int = (9.72*599.72) *(4* h_par_int(i) + h_teto_int(i) +
h_chao_int(i))*( T_ar_int(i-1) - T_par_int(i-1)); %Taxa de calor
trocado por conveccao pela parede interna com o ar interno
841 q_cond = S*k_la*( T_par_ext(i-1) - T_par_int(i-1)); %Taxa de
calor trocado por conducao pela parede interna com a parede externa
842 q_conv_par_ext = (10*600) *(4* h_par_ext(i) + h_teto_ext(i) +
s_chao_ext(i))*(T_amb - T_par_ext(i-1)); %Taxa de calor trocado por
conveccao pela parede externa com o ambiente externo
843 q_rad_par_ext = E*sigma *(2*600*10+2*35*10+600*35) *(T_amb ^4 -
T_par_ext(i-1)^4); %Taxa de calor trocado por radiacao pela parede
externa com o ambiente externo
844
845 m_neve = m_neve + Q*d_agua*dt;
846 q_solid = Q*d_agua*dt*( cp_agua *25+ h_solid);
53
847 q_neve = Q*d_agua*dt*cp_neve *20;
848
849 dT_ar_int = (-q_solid+q_neve+q_lamp + q_evap - q_conv_cubo -
q_conv_par_int)/(m_ar*cp_ar); %Taxa de variacao da temperatura do ar
interno
850 dT_cubo = (q_conv_cubo + q_rad_cubo)/(N*m_cubo*cp_cubo); %Taxa
de variacao da temperatura do cubo
851 dT_par_int = (q_solid -q_neve+q_conv_par_int - q_rad_cubo +
q_cond)/C_int; %Taxa de variacao da temperatura das paredes internas
852 dT_par_ext = (q_rad_par_ext + q_conv_par_ext - q_cond)/C_ext; %
Taxa de variacao da temperatura das paredes externas
853
854 T_ar_int(i) = T_ar_int(i-1) + dt*dT_ar_int; %Temperatura do ar
interno
855 T_cubo(i) = T_cubo(i-1) + dt*dT_cubo; %Temperatura do cubo
856 T_par_int(i) = T_par_int(i-1) + dt*dT_par_int; %Temperatura da
parede interna
857 T_par_ext(i) = T_par_ext(i-1) + dt*dT_par_ext; %Temperatura da
parede externa
858 end
859 end
860
861
862 # Definindo os parametros do controlador
863 setpoint = -5
864 histerese = 2
865
866 # Definindo a lista de tempos (horas)
867 tempos = []
868 hora_atual = datetime.datetime (2023 , 7, 10, 9, 0) # Inicio as 9hrs
869 hora_final = datetime.datetime (2023 , 7, 10, 21, 0) # Fim as 21hrs
870 hora_incremento = datetime.timedelta(minutes =10) # Incremento de 10
minutos
871 while hora_atual <= hora_final:
872 tempos.append(hora_atual)
873 hora_atual += hora_incremento
874
875 # Definindo a lista de temperaturas
876 temperaturas = []
54
877 for tempo in tempos:
878 temperatura = -7 # Temperatura inicial (desligada)
879 if tempo >= datetime.datetime (2023, 7, 10, 12, 0) and tempo <
datetime.datetime (2023 , 7, 10, 13, 0):
880 temperatura = -3 # Temperatura durante a operacao dos canhoes
de neve
881 temperaturas.append(temperatura)
882
883 # Plotando o grafico
884 plt.plot(tempos , temperaturas)
885 plt.xlabel('Tempo ')
886 plt.ylabel('Temperatura (C)')
887 plt.title('Temperaturas ao longo do tempo')
888 plt.xticks(rotation =45)
889 plt.grid(True)
890
891 # Plotando a marcacao para o periodo de operacao dos canhoes de neve
892 plt.axvspan(datetime.datetime (2023 , 7, 10, 12, 0), datetime.datetime
(2023 , 7, 10, 13, 0), facecolor='blue', alpha =0.3)
893
894 plt.show()
895
896 #################################################################
897
898 %QUESTAO 7
899
900 %%%%%% DADOS
901 g = 9.81; %Aceleracao da gravidade
902 sigma = 5.67e-8; %Constante de Stefan -Boltzmann
903 T_amb = 23 + 273.15; %Temperatura doambiente
904 d_ar = 1.17; %Densidade do ar
905 cp_ar = 1003.805; %Calor especifico do ar
906 k_ar = 0.02; %Condutividade termica do ar
907 visc_din = 16e-6; %Viscosidade dinamica do ar
908 visc_cin = visc_din/d_ar; %Viscosidade cinematica do ar
909 Pr = cp_ar*visc_din/k_ar; %numero de Prandtl
910 alfa = k_ar/(d_ar*cp_ar); %Difusividade termica do ar
911 S = 2484.76; %Fator de forma da pista
912 E = 0.90; %Emissividade das superficies
55
913 k_la = 0.033; %Condutividade termica da la de rocha
914 d_la = 32; %Densidade da la de rocha
915 cp_la = 840; %Calor especifico da la de rocha
916 k_ci = 0.35; %Densidade da placa cimenticia
917 d_ci = 1700; %Densidade da placa cimenticia
918 cp_ci = 800; %Calor especifico da placa cimenticia
919 V_int = 599.72*34.72*9.72; %Volume interno da camara
920 V_ext = 600*35*10; %Volume externo da camara
921 V_ci_int = 2*599.72*34.72*0.02+2*34.72*9.72*0.02+2*599.72*34.72*0.02; %
Volume do
922 revestimento interno de aco
923 V_ci_ext = 2*600*10*0.02+2*35*10*0.02+2*600*35*0.02; %Volume do
revestimento
924 externo de placa cimenticia
925 m_la = (V_ext - V_int)*d_la; %Massa da la de rocha
926 m_ci_int = V_ci_int*d_ci; %Massa do revestimento interno de placa
cimenticia
927 m_ci_ext = V_ci_ext*d_ci; %Massa do revestimento externo de placa
cimenticia
928 C_int = m_ci_int*cp_ci + cp_la*m_la /2; %Capacidade termica da parede
interna
929 C_ext = m_ci_ext*cp_ci + cp_la*m_la /2; %Capacidade termica da parede
externa
930 m_evap = 70;
931 cp_evap = 1000;
932 N = 1; %numero de CUBOS
933 A_rad_cubo = 2*2*5; %Area de radiacao dO CUBO DE GELO
934 V_cubo = 8; %Volume do cubo de gelo
935 A_cubo = 2*2*6; %Area dO CUBO DE GELO
936 m_cubo = 7360; %Massa dO CUBO DE GELO
937 cp_cubo = 2100; %Calor especifico dO CUBO DE GELO
938 D_cubo = 1.24; %Diametro da aproximacao dO CUBO para uma esfera
939 m_ar = (V_int - V_cubo)*d_ar; %Massa de ar dentro da camara
940 UA_ligado = 10000; %Coeficiente global do evaporador ligado
941 UA_desl = 200;
942 T_desl = 273.15 -7; %Temp de desligadmento evaporador ( -7 C )
943 T_liga = 273.15 -3; %Temp de religamento do evaporador ( -3 C )
944 % T_evap (1) = 273.15 -20; %Temperatura do evaporador ligado
945 E_gasta_J = 0; %Energia gasta , em Joules
56
946 ra_cubo = g*( D_cubo ^3)/(alfa*visc_cin); %numero de Rayleigh para cubo
sem temperatura
947 ra_par_int = g*(9.72^3) /(alfa*visc_cin); %numero de Rayleigh para parede
interna sem
948 temperatura
949 ra_hor_int = g*(16.41^3) /(alfa*visc_cin); %numero de Rayleigh para
superficies horizontais
950 internas sem temperaturas
951 ra_par_ext = g*(10^3) /(alfa*visc_cin); %numero de Rayleigh para parede
externa sem
952 temperatura
953 ra_hor_ext = g*(16.54^3) /(alfa*visc_cin); %numero de Rayleigh para
superficies horizontais
954 externas sem temperatura
955 tal = 4*60; %4 min
956 q_resist = 5000;
957 %%%%% VARIVEIS DE CONTROLE
958 Evap =1;
959 Ciclos =-1;
960 Degelo =0;
961 int_deg =60*60;
962 t_deg =0;
963 dur_deg =20*60;
964 t_dur_deg =0;
965 %%%%% TEMPO
966 dt=0.1;
967 t_fin =48000;
968 t=0:dt:t_fin;
969 %%%%%% VETORES
970 T_ar_int = zeros (1,(1 + t_fin/dt)); %Temperatura do ar interno
971 T_cubo = zeros (1,(1 + t_fin/dt)); %Temperatura do cubo
972 T_par_int = zeros (1,(1 + t_fin/dt)); %Temperatura da parede interna
973 T_par_ext = zeros (1,(1 + t_fin/dt)); %Temperatura da parede externa
974 h_cubo = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao ddo cubo
975 h_par_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao da parede
interna
976 h_par_ext = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao da parede
externa
57
977 h_chao_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do chao
interno
978 s_chao_ext = zeros (1,(1 + t_fin/dt)); %fator de forma do chao externo
979 h_teto_int = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do teto
interno
980 h_teto_ext = zeros (1,(1 + t_fin/dt)); %Coeficiente de conveccao do teto
externo
981 %%%%% ITERACOES
982 for i = 1:(1 + t_fin/dt)
983 if i ==1 %Condicoes iniciais
984 T_ar_int(i) = -5+273.15;
985 T_cubo(i) = 273.15;
986 T_par_int(i) = T_amb;
987 T_par_ext(i) = T_amb;
988 T_s(i) = T_amb;
989 T_evap(i) = 273.15 -20; %Temperatura do evaporador ligado
990 q_evap = 0;
991 q_lamp = 0;
992 else
993 if (t_dur_deg >= dur_deg)
994 Degelo = 0;
995 t_dur_deg = 0;
996 end
997 if (t_deg <int_deg) && (Degelo ==0)
998 if (Evap ==1) && (T_s(i-1) < T_desl)
999 Evap =0;
1000 Ciclos=Ciclos +1;
1001 if Ciclos ==4
1002 t_inicio = dt*(i-2);
1003 end
1004 if Ciclos == 5
1005 t_final = dt*(i-2);
1006 end
1007 end
1008 if(Evap ==0) && (T_s(i-1) > T_liga)
1009 Evap =1;
1010 T_evap(i-1)=T_evap (1);
1011 end
1012 end
58
1013 if(t_deg >= int_deg)
1014 Degelo =1;
1015 Evap =0;
1016 t_deg =0;
1017 end
1018 Tf_cubo = (T_cubo(i-1) + T_ar_int(i-1))/2; %Temperatura de filme do cubo
1019 Tf_par_int = (T_par_int(i-1) + T_ar_int(i-1))/2; %Temperatura de filme
das paredes
1020 internas
1021 Tf_par_ext = (T_par_ext(i-1) + T_amb)/2; %Temperatura de filme das
paredes externas
1022 B_cubo = 1/ Tf_cubo; %Beta do cubo
1023 B_par_int = 1/ Tf_par_int; %Beta das paredes internas
1024 B_par_ext = 1/ Tf_par_ext; %Beta das paredes internas
1025 Ra_cubo = ra_cubo*B_cubo*abs(T_cubo(i-1) - T_ar_int(i-1)); %numero de
Rayleigh
1026 do cubo
1027 Ra_par_int = ra_par_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(i-1)); %
numero de
1028 Rayleigh das paredes internas
1029 Ra_hor_int = ra_hor_int*B_par_int*abs(T_par_int(i-1) - T_ar_int(i-1)); %
numero de
1030 Rayleigh das superficies horizontais internas
1031 Ra_par_ext = ra_par_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb); %numero
de
1032 Rayleigh das paredes externas
1033 Ra_hor_ext = ra_hor_ext*B_par_ext*abs(T_par_ext(i-1) - T_amb); %numero
de
1034 Rayleigh das superficies horizontais externas
1035 Nu_cubo = 2 + 0.589*( Ra_cubo ^(1/4))/((1 + (0.469/ Pr)^(9/16))^(4/9)); %
numero de
1036 Nusselt do cubo
1037 Nu_par_int = (0.825 + 0.387*( Ra_par_int ^(1/6))/((1+(0.492/ Pr)^(9/16))
^(8/27)))^2;
1038 %numero de Nusselt das paredes internas
1039 Nu_par_ext = (0.825 + 0.387*( Ra_par_ext ^(1/6))/((1+(0.492/ Pr)^(9/16))
^(8/27)))^2;
1040 %numero de Nusselt das paredes externas
59
1041 Nu_teto_int = 0.27*( Ra_hor_int ^(1/4)); %numero de Nusselt do teto
interno
1042 Nu_teto_ext = 0.27*( Ra_hor_ext ^(1/4)); %numero de Nusselt do teto
externo
1043 if Ra_hor_int < 10^7
1044 Nu_chao_int = 0.54*( Ra_hor_int ^(1/4)); %numero de Nusselt do chao
interno para
1045 escoamento laminar
1046 else
1047 Nu_chao_int = 0.15*( Ra_hor_int ^(1/3)); %numero de Nusselt do chao
interno para
1048 escoamento turbulento
1049 end
1050 if Ra_hor_ext < 10^7
1051 Nu_chao_ext = 0.54*( Ra_hor_ext ^(1/4));
1052 else
1053 Nu_chao_ext = 0.15*( Ra_hor_ext ^(1/4));
1054 end
1055 h_cubo(i) = Nu_cubo*k_ar/D_cubo; %Coeficiente de conveccao do cubo
1056 h_par_int(i) = Nu_par_int*k_ar /9.72; %COeficiente de conveccao das
paredes internas
1057 h_par_ext(i) = Nu_par_ext*k_ar /10; %Coeficiente de conveccao das paredes
externas
1058 h_teto_int(i) = Nu_teto_int*k_ar /16.41; %Coeficiente de conveccao do
teto interno
1059 h_teto_ext(i) = Nu_teto_ext*k_ar /16.54; %Coeficiente de conveccao do
teto externo
1060 h_chao_int(i) = Nu_chao_int*k_ar /16.41; %Coeficiente de conveccao do
chao interno
1061 s_chao_ext(i) = Nu_chao_ext*k_ar /16.54;
1062 if Evap ==1
1063 q_evap=UA_ligado *( T_evap(i-1)-T_ar_int(i-1));
1064 elseif Evap ==0
1065 q_evap=UA_desl *( T_evap(i-1)-T_ar_int(i-1));
1066 end
1067 q_conv_cubo = h_cubo(i)*N*A_cubo *( T_ar_int(i-1) - T_cubo(i-1)); %Taxa de
calor
1068 trocado por conveccao pelo cubo com o ar interno
60
1069 q_rad_cubo = E*sigma*A_rad_cubo *( T_par_int(i-1)^4 - T_cubo(i-1)^4); %
Taxa de
1070 calor trocado por radiacao pelo cubo com as paredes internas
1071 q_conv_par_int = (9.72*599.72) *(4* h_par_int(i) + h_teto_int(i) +
1072 h_chao_int(i))*( T_ar_int(i-1) - T_par_int(i-1)); %Taxa