UNIVERSIDADE ESTADUAL DE CAMPINAS FACULDADE DE ENGENHARIA
MECNICA TRABALHO DE GRADUAO II - EM919 Uma Anlise Terica do Motor
Stirling Tipo Beta Orientador: Pablo Siqueira Meirelles Autor:
Felipe Santin Furlan RA: 043429 Campinas Junho de 2011 ii
Agradecimentos
Atodasaspessoasqueapoiaram,emespecialaDouglasGeraldiquesugeriua
idia do motor Sitlring e ao professor Pablo Siqueira Meirelles que
incentivou e orientou o assunto. A motivao deste trabalho surgiu de
participao de uma competio de projetos,
naqualcriamosumprottipodomotorStirling.Emboraatomomentoelenotenha
funcionado,essetrabalhoservirdeapoioparaevidenciarosfenmenosetambm
para ajudar futuramente na construo de novos prottipos. iii Resumo
Esse trabalho elabora uma equao de movimento do motor Stirling tipo
beta com
nfasenomecanismo,termodinmicaetransfernciadecalor.Nasdimenses
utilizadas o motor no chega a atingir 1watts, operando com
temperaturas na superfcie
de50Cnapartefriaede400Cnapartequente.Vamosanalisaravariaode
temperatura do fluido em cada volume de controle do motor, a
influncia da quantidade de massa na gerao potncia e a influencia do
torque do desempenho do motor. iv Sumrio
Introduo...............................................................................................................
6 Reviso
Bibliogrfica...............................................................................................
7
Metodologia............................................................................................................11
Mecanismos
.......................................................................................................14
Mecanismo de
Potncia..................................................................................15
Mecanismo de
Sincronia.................................................................................23
Mecanismo de
Deslocamento.........................................................................27
Foras
Generalizadas.....................................................................................32
Energia Potencial
............................................................................................33
Termodinmica e Transferncia de Calor
..........................................................34 Calor
...............................................................................................................38
Trabalho..........................................................................................................38
Presso...........................................................................................................39
Anlise de
Desempenho.....................................................................................40
Resultados
.............................................................................................................41
Concluso..............................................................................................................54
Referncias
Bibliogrficas......................................................................................55
v ANEXO I Tabela das Propriedades do AR como gs
ideal..........................56 ANEXO II Clculo dos Coeficientes
de Conveco .....................................57 ANEXO III Inrcia
dos Elementos
................................................................58
ANEXO IV Mecanismo de Sincronia: calculo dos
ngulos...........................59 ANEXO V Programa em
Matlab...................................................................61
6 Introduo
Em1816opastorescocsRobertStirlingpatenteouomotorStirling.Entretanto
at hoje no tem uma utilizao massiva. Sua idia simples e seu
potencial terico de alta eficincia estimulam desde ento centros de
pesquisas, universidades e empresas privadas no seu desenvolvimento
e aprimoramento. Ele apresenta caractersticas interessantes como
flexibilidade. Por exemplo, com relao entrada de energia, ele pode
operar com qualquer fonte quente; com relao
aociclotermodinmico,podemosterummotorquandocriamosdiferenasde
temperatura, ou um refrigerador ou aquecedor quando fornecemos
potncia.Outracaractersticainteressantequeseuaquecimentoexterno,portanto,se
operado com combustveis, podemos controlar a combusto para realizar
umaqueima perfeita e assim, reduzir a emisso de elementos
prejudiciais ao meio ambiente. Nas ultimas dcadas ele foi mais
utilizado como refrigerador de cmara criognica
doquecomomotor.Entretantocomforteapeloambientaleodesenvolvimento
tecnolgico,algumasempresasestosurgindonomercadocomercializandoomotor
Stirling movido energia solar acoplado a uma gerador de energia
eltrica.[4].
EssetrabalhoirrealizarumaanlisetericadofuncionamentodomotorStirling
com nfase em seus mecanismos e seu ciclo termodinmico. 7 Reviso
Bibliogrfica Atualmente temos inmeros sistemas que realizam
trabalho mecnico, na figura 1
abaixopodemosclassificaromotorStirlingdeacordocomalgumasespecificaese
diferencia-lo de outros motores. Figura 1 - Classificao do Motor
Stirling O motor Stirling uma mquina trmica no qual fluidos gasosos
realizam trabalho
emumclicofechadoporaquecimentoexterno.Temosdiversosmodelos,entretanto
uma classificao usual dividi-los em tipo Alfa, Beta e Gama, como
exemplificado na figura 2. 8 Figura 2 - Modelos de motor Stirling
Nesse trabalho iremos analisar o motor tipo Beta, entretanto todos
eles realizam o
mesmociclotermodinmicoqueapresenta4faseseexecutadoem2temposdo
pisto:compressoisotrmica(temperaturaconstante),aquecimentoisovolumtrico
(volumeconstante),expansoisotrmica(temperaturaconstante)eresfriamento
isovolumtrico (volume constante), como demonstrado na figura 3. 9
Figura 3 - Ciclo Termodinmico Stirling (baseado em DARLINGTON[4])
Ofuncionamentodesteciclotemalgunsdetalhesimportantesquesero
exemplificados na figura 4. 10 Figura 4 - Descrio do funcionamento
do Motor Stirling 11 Metodologia
OfuncionamentodomotorStirlingconseqnciadapressoresultantedo
aquecimentoedoresfriamentodofluidodetrabalho.Apressoatuanopistode
potncia que movimenta o
mecanismo.Paradeterminaraequaodemovimentodosistemavamosutilizarasegunda
forma geral de Lagrange que leva em considerao a energia potencial
do sistema.ncjj jQqLqLdtd=&(1.0) Otermo ncjQ
compostoporumvetorcolunaquerepresentaasforas
generalizadasdevidosforasnoconservativasaplicadasnosistema.Avarivel
jqrepresentatodasascoordenadasindependentesdosistema.SendoLachamada
funoLagrangeanaequetemporfinalidaderepresentaradiferenaentreenergia
cintica e energia potencial do sistema analisado, temos: Ep T L =
(1.1) Sendo que T representa a energia cintica e Ep representa a
energia potencial do sistema. Aps substituir a Equao 1.1 em 1.0
obtm-se: ncjj j j jQqEpqTqEpdtdqTdtd=+& & (1.3) A energia
potencial de um sistema obtida apenas em funo da coordenada de
posio do mesmo. A partir disso, tem-se: 12 0 =jqEp&(1.4)
Comisso,pode-sesimplificaraequao1.3,vistoqueaenergiapotencialno
depende da velocidade do sistema.
ConformeDoughty[5],aenergiacinticadeumsistemapodeserrepresentada
como: ] [ ] [ ] [ ] [ ] [ 5 . 0jT Tjq Kc M Kc q T & & =
(1.5) No qual [jq& ] o vetor velocidades generalizadas, [Kc] a
matriz de coeficientes de velocidades e [M] a matriz de inrcia do
sistema.Temos ento: ((((
=jn jnqqN q qqT&&& & ... ] [ ] ... [11No qual:j
n 1 (1.6) Sendo] [nN : ]) [ ] [ ] [ ] [ ] [ ] ([ 5 , 0 ] [nT Tn nL
M Kc Kc M L N + = (1.7) No qual:] [ ] [ KcqLnn=(1.8)
Apscalcularaderivadadaenergiacinticaemrelao jq& ,derivou-sea
expressoresultanteemrelaoaotempo.Assimdepoisdealgumasoperaes,
obteve-se: 13 ((((
+ + +((((
=jj jTjqqN q N qn qqKc M KcqTdtd&&& &&
&& &&... ]) [ 2 .... ] [ 2 ( ... ] [ ] [ ] [11
11(1.9)
Assumindoquenotemosfolganoseixos,nomotorStirlingtemosapenasuma
varivel independente que chamamos de 1 , dessa forma, a equao que
representa o movimento do sistema : ncjTQEpN N Kc M Kc =+ + 11 1 1
1 1 1 1) ] [ ( ) ] [ 2 ] [ ] [ ] [ ] ([ & & & &
& &(1.10) E podemos simplificar para: ncjTQEpN Kc M Kc =+ +
11 1 1 1] [ ] [ ] [ ] [ ] [ & & & &(1.11)
Paraencontraraconstante[Kc],otermo] [1 N eEpprecisamosresolvera
equaes de loop dos mecanismos do sistema. No termo ncjQ temos a
presso do fluido
detrabalhoquesercalculadopelosfundamentosdetransfernciadecalore
termodinmica. 14 Mecanismos O motor Sitrling tipo beta tem diversas
peas em movimento. O conjunto de peas exemplificadas em azul na
figura 5 ser chamada mecanismo de potncia, pois recebe
apressodofluidodetrabalho;oconjuntodepeasemverdetemonomede
mecanismo de sincronia, pois mantm entre os pistes o ngulo de fase
determinado; o conjunto de peas em vermelho recebe o nome de
mecanismo de deslocamento, pois responsvel pelo deslocamento do
fluido de trabalho entre a cmara fria e a cmara quente. Figura 5 -
Mecanismo: em azul de potencia, em verde de sincronia, em vermelho
de deslocamento 15 Mecanismo de Potncia Figura 6 - Mecanismo de
Potencia Percorrendo o loop em x: 0 4 1 5 ) cos( 4 ) cos( 14 1= + +
M M M Xp B B (2.0) Percorrendo o loop em y, considerando M3 igual
M2: 0 ) ( 4 ) ( 14 1= + sen B sen B(2.1) Com as duas equaes podemos
isolar 4 eXp , assim temos: ||
\| =4) ( * 11 14Bsen Bsen(2.2) 4 1 54) ( * 1cos 4 ) cos( 14 1 5
) cos( 4 ) cos( 11 114 1M M MBsen Bsen B B XpM M M B B Xp +|||
\|||
\| + = + + = (2.3) 16 Derivando no tempo 2.0 e 2.1 temos: 0 ) (
4 ) ( 14 4 1 1= p X sen B sen B& & & (2.4) 0 ) cos( 4 )
cos( 14 4 1 1= + & &B B(2.5) Isolando 4& temos: 1414)
cos( 4) cos( 1& & =BB(2.6) Temos que o coeficiente de
velocidade de 4&:
) cos( 4) cos( 1414 =BBKc (2.7) Isolandop X& temos: 4 4 1 1)
( 4 ) ( 1 & & & = sen B sen B p X (2.8) Logo: 4 4 1 1)
( 4 ) ( 1 Kc sen B sen B Kcp = & (2.9) Derivando no tempo 2.4 e
2.5 temos: 0 ) cos( 4 ) ( 4 ) cos( 1 ) ( 124 4 4 421 1 1 1= p X B
sen B B sen B& & & & & & & & (2.10)
0 ) sin( 4 ) cos( 4 ) sin( 1 ) cos( 124 4 4 421 1 1 1= + &
& & & &B B B B(2.11) Isolando 4& & temos:
17 1 4 1 4 4424 4 4 421 1 1 14) cos( 4) sin( 4 ) cos( 4 ) sin( 1 )
cos( 1 & & & & && & & & &
&& &L KcBB B B B+ = + + + = (2.12) No qual 4L o
coeficiente de acelerao de 4& &. Isolando pX&
&temos: 24 421 1 4 4 1 1) cos( 4 ) cos( 1 ) ( 4 ) ( 1 &
& & & & & & & = B B sen B sen B p X
(2.13) ( )21 12124 421 1 1 4 4 1 1) cos( 4 ) cos( 1 ) ( 4 ) ( 1
& & & & && & & & & &
& & + = =p pL Kc p XKc B B Kc sen B sen B p X No qual pL o
coeficiente de acelerao dep X& &. Centro de Massa das peas
do mecanismo de potencia
SendoocentrodemassadeB1opontoPb1comcoordenadalocalemuPb1e
vPb1,vamostransformaremcoordenadoglobaleencontrarasconstantesde
velocidade e acelerao para cada pea do mecanismo. Figura 7 - Centro
de Massa de B1 Em x temos: 18 1 121 1 1 121 1 11 1 1 1 11 1 1) cos(
1 ) ( 1 ) ( 1 ) cos( 1) cos( 1 ) ( 1) ( 1 ) cos( 1 & &
& & & & & && & & + = = =vPb sen
vPb sen uPb uPb XvPb sen uPb Xsen vPb uPb XPbPbPb Dessa forma
encontramos: ) ( 1 ) cos( 1) cos( 1 ) ( 11 1 11 1 1 sen vPb uPb
LvPb sen uPb KcXPbXPb + = =(2.16) Em y temos: 1 1 1 1 11 1 1) ( 1 )
cos( 1) cos( 1 ) ( 1 & & & = + =sen vPb uPb YvPb sen
uPb YPbPb(2.17) Dessa forma encontramos: ) ( 1 ) cos( 11 1 1 sen
vPb uPb KcYPb = (2.18) Como ( )111 dKc dLYPbYPb= temos: ) cos( 1 )
( 11 1 1 = vPb sen uPb LYPb(2.19) Da mesma forma vamos encontra as
constantes de velocidade e acelerao para a pea B4. 19 Figura 8 -
Centro de massa de B4 Em x temos: 4 424 44 424 4 1 121 1 44 4 4 4 1
1 44 4 1 4) cos( 4 ) ( 4) ( 4 ) cos( 4 ) ( 1 ) cos( 1) cos( 4 ) ( 4
) ( 1) ( 4 ) cos( 4 ) cos( 1 & & && & &
& & & & && & & & ++ = = + =vPb
sen vPbsen uPb uPb sen B B XvPb sen uPb sen B Xsen vPb uPb B
XPbPbPb Logo as constantes so: ( )( ) ( )24 4 4 1 44 4 4 1 4) cos(
4 ) ( 4 ) cos( 1) cos( 4 ) ( 4 ) ( 1 Kc uPb sen vPb B LKc vPb sen
uPb sen B KcXPBXPB + = + =(2.21) Em y temos: 4 4 4 4 1 1 44 4 1 4)
( 4 ) cos( 4 ) cos( 1) cos( 4 ) ( 4 ) sin( 1 & & &
& + = + =sen vPb uPb B YvPb sen uPb B YPbPb(2.22) 20 Logo as
constantes so: ( )( ) ( )24 4 4 1 44 4 4 1 4) cos( 4 ) ( 4 ) ( 1) (
4 ) cos( 4 ) cos( 1 Kc vPb sen uPb sen B LKc sen vPb uPb B KcYPBYPB
+ = + =(2.23) A seguir repetimos o procedimento para encontrar as
constantes para o Pisto de Potncia. Figura 9 - Centro de Massa do
Pisto de Potncia Em x: 4 424 4 1 121 14 4 1 14 1) ( 4 ) cos( 4 ) (
1 ) cos( 1) ( 4 ) ( 1) cos( 4 ) cos( 1 & & & &
& && & & = =+ + =sen B B sen B B Xsen B sen B
XuPxp B B XPXPPXPPXP (2.24) Assim obtemos: ( )24 4 14 1) cos( 4 )
cos( 1) ( 4 ) ( 1 Kc B B Lsen B sen B KcXPXPXPXP = =(2.25) Como o
pisto no se move em y, as constantes so igual a zero, por isso no
so necessria. 21
NessemecanismotemosacopladoaB1ocontrapesoeodiscodeinrciaque ameniza
a variao de velocidade do motor. Iremos fazer uma analise dessas
peas de forma separada de B1. .Figura 10 - Disco de Inrcia, contra
peso e biela 1
ApesardaBilea1serdimensionadaportodoocontornoemvermelhonafigura
10,iremosdividi-laemtrspartes.AbielaB1emazulclaro,ocontrapesoemazul
escuroeumapartedesprezadaembranco.Emverdeodiscodeinrcia.Issofeito
parasimplificaroclculodeinrcia,bemcomoanalisarainfluenciadecadapeano
mecanismo. AsconstantesderotaoetranslaodapeaContraPeso(CP)soas
seguintes: 22 ; 0; 1) 1 cos( ) 1 sin() 1 ( ) 1 cos() 1 ( ) 1 cos()
1 cos( ) 1 (11== = + = = =CPCPYCPXCPYCPXCPLKcbeta vCP beta uCP
Lbeta sen vCP beta uCP Lbeta sen vCP beta uCP Kcbeta vCP beta sen
uCP Kc (2.26) E para o disco de inrcia temos apenas as constantes
de rotao: 01==DIDILKc(2.27) 23 Mecanismo de Sincronia
NomecanismorealB3eB5soumapeanica,entretantoparaanlisedo mecanismo
dividimos essa pea em duas, como pode ser visto na figura 10.
Figura 11 - Diviso de uma pea em duas O mesmo procedimento ser
feito para o mecanismo de sincronia para encontrar as constantes de
velocidade e acelerao de rotao e de translao. Figura 12 - Mecanismo
de Sincronia 24 Em x2: 0 4 ) cos( 1 ) cos( 2 ) cos( 31 2 3= + + M B
B B (2.28) Em y2: 0 2 ) ( 1 ) ( 2 ) ( 31 2 3= + M sen B sen B sen B
(2.29) Depois de alguns passos, demonstrado no ANEXO IV chegamosa
uma equao de 2 grau da forma: 0 ) ( cos ) ( cos2222= + + c b a
(2.30) Na qual: ( ) ( )( ) ( )( ) ( )2 3 22 3) ( cos) ( 3 4 ) cos(
12 ) cos( 3 4 ) cos( 1 22 ) cos( 3 ) ( 32 2 2 212 2112 2B BC D B
Bww zz sen B M B cB z B M B bB z B z sen B a + + == =+ =+ + =(2.31)
Por conseguinte:z =2 3 (2.32) Derivando temos: 0 cos 1 cos 2 cos 30
1 2 31 1 2 2 3 31 1 2 2 3 3= + = + & & && &
&B B Bsen B sen B sen B(2.33) Encontramos as constantes de
velocidade: 25 ) cos( 3) cos( 1 ) cos( 2) ( cos ) () tan( cos2131 2
232 2 31 3 12 + =((
=BB Kc BKcsen tgsenBBKc (2.34) E encontramos as constantes de
acelerao: ( ) ( ) ( )( )( ) ( ) ( )) ( 31 2 32cos 1 cos 2 cos 33122
223 332122 223 32 sen Bsen B Kc sen B Kc sen BLsen BB Kc B Kc BL +
= + =(2.35) Centro de Massa do Mecanismo de Sincronia
Omecanismodesincroniaresponsvelpeloongulodefaseentreo
mecanismodepotnciaededescolamento.Nonossocasoemestudo,ongulode fase
constante. Figura 13 - Centro de massa do mecanismo de sincronia 26
) ( 3 ) cos( 2 ) ( 2) cos( 3 ) ( 2 ) cos( 2) cos( 3 ) ( 3 ) 180
cos( 3 ) 180 ( 3) ( 3 ) cos( 3 ) 180 ( 3 ) 180 cos( 33 2 2 23 2 2
23 3 3 3 33 3 3 3 3 sen B vPb sen uPb YB sen vPb uPb XvPb sen uPb
vPb sen uPb Ysen vPb uPb sen vPb uPb XPBPBPBPB + + = + = = + = + =
=(2.37) Derivandonotempoerealizandoalgumasmanipulaesencontramosas
constates de velocidade: ( )( )( )( )3 3 2 2 2 23 3 2 2 2 23 3 3 33
3 3 3) cos( 3 ) ( 2 ) ( 2) ( 3 ) cos( 2 ) ( 2) ( 3 ) cos( 3) cos( 3
) ( 3 Kc B Kc sen vPb sen uPb KcKc sen B Kc vPb sen uPb KcKc sen
vPb uPb KcKc vPb sen uPb KcYPBXPBYPBXPB + = = + = + =(2.38)
Derivando novamente encontramos as constantes de acelerao: ( )( )(
)( )23 322 2 2 223 322 2 2 223 3 3 323 3 3 3) ( 3 ) cos( 2 ) cos(
2) cos( 3 ) ( 2 ) cos( 2) cos( 3 ) ( 3) ( 3 ) cos( 3 Kc sen B Kc
vPb uPb LKc B Kc sen vPb uPb LKc vPb sen uPb LKc sen vPb uPb
LYPBXPBYPBXPB = + = + = =(2.39) 27 Mecanismo de Deslocamento Figura
14 - Mecanismo de Deslocamento Nessemecanismotemosoelofixoentre 5
(beta5)e 3 (beta3),nocaso A (betaA).A + =3 5(2.40) Fazendo o loop
de posio em x e em y, respectivamente: 0 2 6 50 1 7 cos 6 cos 56 56
5= + + = + + M sen B sen BM Xd B B B (2.41) Derivando temos: 0 cos
6 cos 50 6 56 6 5 56 6 5 53 5= + = + = & && &
&& &B Bd X sen B sen B (2.42) 28 As constantes de
velocidade so: 6 6 5 556563 5) ( 6 ) ( 5) cos( 6) cos( 5 Kc sen B
Kc sen B KcKcBBKcKc KcXD = ==(2.43) As constantes de acelerao so: (
) ( )( ) ( )26 625 5626 625 563 5) cos( 6 ) cos( 5cos 6) ( 6 ) ( 5
Kc B Kc B LBKc sen B Kc sen BLL LXD = + ==(2.44) Centro de Massa
Mecanismo de Deslocamento Figura 15 - Centro de Massa mecanismo de
Deslocamento 29 27 ) cos( 6 ) cos( 5) cos( 6 ) ( 6 ) ( 5) ( 6 )
cos( 6 ) cos( 5) cos( 5 ) ( 5) ( 5 ) cos( 56 56 6 5 66 6 5 65 5 55
5 5M YuPxd B B B XvPb sen uPb sen B Ysen vPb uPb B XvPB sen uPb
Ysen vPb uPb XPDPDPBPBPBPB=+ + + = + + = + = + = = (2.45) A partir
dessas equaes nos derivamos no tempo e encontramos as constates de
velocidade. 0) cos( 6 ) ( 5) ( 6 ) cos( 6 ) cos( 5) cos( 6 ) ( 6 )
( 5) ( 5 ) cos( 5) cos( 5 ) ( 56 6 5 56 6 6 6 5 5 66 6 6 6 5 5 65 5
5 5 55 5 5 5 5= + = + = = = =YPDXPDYPBXPBYPBXPBKcKc B Kc sen B KcKc
sen vPb Kc uPb Kc B KcKc vPb Kc sen uPb Kc sen B KcKc sen vPB Kc
uPb KcKc vPb Kc sen uPb Kc (2.46)
Derivandoasconstantesdevelocidadenotempoencontramosasseguintes
constantes de acelerao: ( ) ( )( ) ( )( ) ( ) ( )( ) ( ) ( )( ) (
)0) ( 6 ) cos( 5) cos( 6 ) ( 6 ) ( 5) ( 6 ) cos( 6 ) cos( 5) cos( 5
) ( 5) ( 5 ) cos( 526 625 525 626 625 5 626 626 625 5 625 525 5 525
525 5 5= = = = = + =YPDXPDYPBXPBYPBXPBLKc sen B Kc B LKc vPb Kc sen
uPb Kc sen B LKc sen vPb Kc uPb Kc B LKc vPB Kc sen uPb LKc sen vPb
Kc uPb L (2.47) 30 Matrizes Depois de calculado todas a constante
de velocidade e acelerao que envolve o nosso sistema, vamos
organiz-las em forma de matrizpara posteriormenteaplicar na equao
de movimento. Algumas constantes esto duplicadas, como por exemplo,
a do pisto e outros so nulos, assim j sero eliminadas dessa matriz.
1654321665544332211654321665544332211] [111
=(((((((((((((((((((((((((((((((
=((((((((((((((((((((((((((((((((((
KcKcKcKcKcKcKcKcKcKcKcKcKcKcKcKcKcKcKcKcKcKcKcKcKcYXXXYXYXYXYXYXYXDICPYCPXCPXPDXPPYPBXPBYPBXPBYPBXPBYPBXPBYPBXPBYPBXPBDIXPCPCPPDPPPBPBPBPBPBPBPBPBPBPBPBPB&&&&&&&&&&&&&&&&&&&&&&&&
(((((((((((((((((((((((((((((((
=DICPYCPXCPXPDXPPYPBXPBYPBXPBYPBXPBYPBXPBYPBXPBYPBXPBLLLLLLLLLLLLLLLLLLLLLLLLL1654321665544332211
31 Seguindooordenamentodamatrizdevelocidadeeacelerao,temosquea
matriz de inrcia :((((((((((((((((((((((((((((((
=DICPCPCPPDPPBBBBBBBBBBBBMIMMIIIIIIMMMMMMMMMMMMMMM110 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0653321665544332211
NoqualMxxamassadoelementoXXeIxxainrciadoelementoXXem
relaoaoseucentrodegiro.Oselementosforamsimplificadosparaformas
geomtricassimplesparafacilitarocalculodomomentodeinrcia.NoANEXOIII
temos as forma utilizadas para cada elemento. 32 Foras
Generalizadas
Todaforaatuanteemumdeterminadosistemateminfluncianaresposta dinmica
do mesmo. Com isso, necessrio determinar uma nica fora generalizada
Qque,quandoatuandoatravsdeumavariaodacoordenadavirtualq,farum
trabalhovirtualQ.q,igualsomadostrabalhosvirtuaisdasforasatuantes,quese
deslocam atravs dos seus deslocamentos virtuais associados
(Doughty, 1988). A partir disso pode-se escrever:q Q r F Wi i i i =
+ = (2.50) Figura 16 - Foras atuando no sistema A varivel Fp uma
fora aplicada na cabea do pisto. Essa fora representa a
pressodofluido.JavarivelTorque( )umtorqueexternoaplicadonamanivela
querepresentaaresistnciadosistema.Dessaformapodemosescreverotrabalho
virtual do sistema como: + = = + = + =XPPFPKc Fp QQddXFp X Fp W1 1
111 (2.51) 33 Sendo que Q representa a fora generalizada do
sistema.
Aforaestrelacionadacomapresso(p)resultantedofluidoserobtidano
tpicoTERMODINAMICAETRANSFERENCIADECALOR.Assim,comapressop
calculada, a fora que atua no pisto de potencia : 42Dp Fp =(2.52)
No qual D o dimetro do cilindro.
Nocasodotorque,vamossimularqueelevarialinearmentecomrelaoa
velocidade. Dessa forma: 1 & =K (2.53) Energia Potencial
Podemos calcular a energia potencial do sistema como: ( )6 6 5 5 4
4 3 3 2 2 1 1 PB B PB B PB B PB B PB B PB BY M Y M Y M Y M Y M Y M
g Ep + + + + + = (2.54) Assim o termo referente energia potencial
fica: ||||
\| + ++ + + + =6 6 6 5 5 54 4 4 3 3 3 2 2 2 1 1 1...... Kc Kc M
Kc Kc MKc Kc M Kc Kc M Kc Kc M Kc Kc MgdqdEpYPB B YPB BYPB B YPB B
YPB B YPB B 34 Termodinmica e Transferncia de Calor Nesse tpico
iremos abordar os fenmenos envolvidos na variao de presso do
sistema.Hipteses simplificadoras: 1.Massa de fluido constante, isto
, no h vazamento do cilindro para o meio externo e do meio externo
para o cilindro.
2.Temperaturadasuperfcieconstanteemcadavolumedecontrole.Dessa
formatemosatemperaturaTsfnasuperfcie(azul)dovolumedecontrole
dapartefria,Tsqfnasuperfcie(amarela)dovolumedecontroledo
regenerador e Tsq na superfcie (vermelha) do volume de controle da
parte quente. Figura 17 - Distribuio de Temperatura no Cilindro
Como estamos tratando com fluidos gasosos no interior do cilindro,
a equao que governa os parmetros levando em conta a
compressibilidade : T R nV pZ = (3.0) 35
NoqualZofatordecompressibilidade,paPresso,Vovolume,Ra constante
universal dos gases, n o numero de Mols no sistema e T a
temperatura.
Considerandoqueofluidosecomporteconformeomodelodeumgsideal,Z tende
a 1, dessa forma temos: T R nV p = 1 (3.1)
Essahiptesesetornavlidaparaestadosondeapressoppequenaem
relaopressocriticapc(baixapr)e/ouatemperaturaTelevadaemrelao
temperaturacrticaTc(elevadaTr)ofatordecompressibilidadeprximode1,e
podemos admitir com um preciso aceitvel que Z=1. (INCROPERA, pg.
77, 2000) [1]
Apesardeapressoseramesmaemtodososvolumesdecontrole,ogsest
distribudonocilindroemdiferentesconcentraes,poisadensidadeaolongodele
varivel. Entretanto, sabemos pela hiptese 1 que a soma das massas
em cada volume
decontroleresultanamassatotal..Damesmaformasomandoaquantidadedemols
em cada volume de controle o resultado a quantidade de mols total,
assim temos: qf f qn n n n + + = (3.2) Substituindo 3.2 em 3.1,
temos: ( ) T R n n n V pqf f q + + = (3.3)
Cadavolumedecontroleestsujeitoasuatemperaturaevolume,entretantoa
presso sobre o gs igual em todos os volumes de controle. Assim
temos: T RT RpVT RpVT RpVV pffqfqfqq |||
\|++= (3.4) 36 Rearranjando 3.4 e isolando T, temos: VV T T V T
T V T TT T TTqf f q f qf q q qf ff qf q + + =(3.5) No caso T a
temperatura mdia do fluido considerando as temperaturas do fluido
emcadavolumedecontrole.AolongodotempoatemperaturamdiadofluidoTeo
volume total do sistema V se altera, alterando tambm a presso p.
Figura 18 - Variao do volume de controle
Vamosconsiderarqueaspropriedadesdofluidopermaneamconstantesno
intervalodetempodet0att1.Dessaformaobalanodeenergiaemumvolumede
controle transiente se torna: + = ss see e VC VC VC VCh m h m W Q U
t U ) 0 ( ) ((3.6) Na qual U a energia interna, Q a quantidade e
energia transferia por calor, W
aquantidadedeenergiatransferidaportrabalhoeaoutrasparcelasrepresentamas
37 somatrias da energia que entram e que sai do volume de controle,
m a massa e h e a entalpia. Estamos desprezando a energia cintica e
potencial do fluido. Figura 19 Balano de Energia nos volumes de
Controle
Sendom1aquantidadedemassaquesaiouqueentranovolumedecontrole quente
e m2 a quantidade de massa que entra ou sai do volume de controle
frio. Temos quando a massa esta entrando nos volumes de controle: (
)( ) )) ( ( ) ( ) ( 2)) ( ( ) ( ) ( 10 0 10 0 1t T t Vf t Vf mt T t
Vq t Vq mqfqf = =(3.7) E quando esta saindo ( )( ) )) ( ( ) ( ) (
2)) ( ( ) ( ) ( 10 0 10 0 1t T t Vf t Vf mt T t Vq t Vq mfq = =
(3.8) No qual)) ( (0t Tq a densidade do fluido na temperatura Tq no
instante t0. O valor
daquantidadedemassapositivoquandoestentrandonovolumedecontrolee
negativo quando est saindo.Deve-se tambm avaliar corretamente a
entalpia quando a massa entra e quando
amassasai.Nocasodemassasaindodovolumequentetemosqueaenergia 38 ) (
1 Tq h me quando est entrando ) ( 1 Tqf h m . De forma anloga
acontece nos outros volumes de controle. Calor O calor positivo
quando entra e negativo quando sai do volume de controle e se
comporta segundo as equaes: ( )( )( ) ) ( ) () ( ) () ( ) (t T T A
hc t Qt T T A hc t Qt T T A hc t Qqf sqf qf qf qff sf f f fq sq q q
q = = = (3.9) Noqual qhc , fhc e qfhc
sooscoeficientesdeconvecodosrespectivos
volumesdecontrole.Ocalculodoscoeficientesdeconvecoestodemonstradono
ANEXO II. Trabalho
Otrabalhorealizadoquandoofluidoseexpandeousecomprime.Elese expande
na cmara quente quando massa est entrando e se comprime na cmara
fria
quandomassatambmestentrando.Vamosdesconsiderarotrabalhorealizadono
regenerador. Assim temos: ( )( ) ) ( ) ( ( ) ( ) () ( ) ( ) ( ) (0
1 0 10 1 0 1t V t V t p t Wprt V t V t p t Wex = =(3.10)
Otrabalhodeexpanso(Wex)zeroquandomassaestasaindodacmara quente e o
trabalho de compresso (Wpr) tambm zero quando massa est saindo da
cmara fria. Vale salientar que na analise de motor, o trabalho
retira energia do fluido.
EntretantoquandooStirlingatuacomoumrefrigeradore/ouaquecedorotrabalho
fornecido ao fluido. 39 Temperatura em cada volume de controle
Paragsidealtemosqueavariaodaenergiainternadependesomenteda
temperatura, conforme a seguinte equao: ( ) ) ( ) ( )) ( ( ) ( ) (1
2 1 1 2t T t T t T c t U t Uv = (3.11) Logo substituindo 3.6 em
3.11 temos para o volume de controle quente: ( )) ()) ( () ( 1 ) (
) () () ( 1 ) ( ) ( ) ( ) ( ) (111 121 1 1 2 1t Tt T cTx h m t Wex
t Qt TTx h m t Wex t Q t T t T t cqq vqqq v+ + = + = (3.12)
NoqualTxigualaTqquandom1estsaindo,TxigualaTqfquandom1est entrando
no volume de controle quente.Substituindo 3.6 em 3.11 e isolando
T(t2) para o volume de controle frio, temos: ) ()) ( () ( 2 ) ( )
() (111 12t Tt T cTy h m t Wpr t Qt Tff vff+ + = (3.13)
NoqualTyigualaTfquandom2estsaindo,TyigualaTqfquandom1est entrando
no volumede controle frio.Substituindo 3.6 em 3.11 e isolando
Tqf(t2)parao volume de controle do regenerador (Vqf), temos: ) ())
( () ( 2 ) ( 1 ) () (1112t Tt T cTy h m Tx h m t Qt Tqff vfqf+ =
(3.14) Presso
Comotemosastemperaturasdecadavolumedecontrole,podemoscalculara
temperatura mdia no instante t com a equao 3.5. Portanto: 40 ) ()
() (t Vt T R nt p = (3.15) No qual) (t T a temperatura mdia no
instante t e) (t V o volume total no instante t. Anlise de
Desempenho A metodologia ser implantada em cdigo Matlab disponvel
no ANEXO V.
Vamosanalisarodesempenhodosistemacomrelaoapotncia,calculadada
seguinte forma: 6026021RPM v K RPMP = = No qual: |||
\| |||
\|=602) ( ) (1 1i fi ft tt tRPM Sendo 1 v
avelocidademdiaeRPMarotaoporminutomdia. 41 Resultados O Motor
Sitrling tipo Beta que vamos simular tem as seguintes grandezas:
Dimenso das Peas [metros]Posio Centro de MassaDescrioDistncia entre
eixos LarguraAlturaComprimentoDireo U Direo V Massa [Kg]
B10.0110.0030.0100.0500.011/200.3 B20.0340.0030.0100.0470.034/200.3
B30.0430.0030.0100.0550.043/200.3 B40.0780.0030.0100.0950.078/200.3
B50.0440.0030.0100.0550.044/200.3 B60.0120.0030.0100.0250.012/200.3
Contra-Peso0.0030.0070.012-0.0080.0130.1765
DescrioValorDescrioValor M10.024Dim. do Cilindro (Dc)0.035
M20.045Comprimento do Cilindro (Cc)0.168 M3M2Compri. do Deslocador
(Cd)0.095 M40.026Dim. do Deslocador (Dd)0.033 M50.016Dim. Disco de
Inrcia0.092 betaA77 Massa do disco de Inrcia3.0 Temperatura na
Parte Quente (Tq)400C Temperatura na Parte Fria (Tf)50C Temperatura
no Regenerador (Tqf)225C
Acondioinicialparaaequaodiferencialdemovimentoposiodebeta1 igual a
zero e velocidade igual a -5rad/s.
VamosutilizarARcomofluidodetrabalho,portantoiremosutilizartabelasdas
propriedades do AR como gs ideal para os clculos termodinmicos e de
transferncia
decalor.Aquantidadedefluidodetrabalhotemumaimportnciavitalparao
funcionamento do motor Stirling. No grfico 1 abaixo, podemos obter
a quantidade mais
eficientedenumerodemolsque7,51E-4,queparaoARsignifica21,32gramas,
42
(nar=28,97kg/mol).Essasimulaomanteveconstantetodososoutrosparmetros,
como torque, temperatura das partes e tempo de simulao. Potncia x
Mols0,000,100,200,300,400,500,600,00E+00 2,00E-04 4,00E-04 6,00E-04
8,00E-04MolsPotncia [Watts] .n Grfico 1 - Potncia x Nmero de Mols
Podemosobservarquetantoaentradacomoasadadeardocilindropode
prejudicar a desempenho do motor. Operando no ponto de mxima
potncia a sada de ar menos prejudicial que a entrada de ar.
Nogrfico2podemosobservarqueexisteumaconstantedetorquetimo.Para
Kt=0,010temos79RPM,potnciade0,65Wattseumtorquede0,0785N*mem100
segundos de simulao. Para as simulaes seguintes usaremos Kt=10. 43
Tabela 1 - Comportamento do motor em funo de Kt RPM e Potncia x
Torque0,000,200,400,600,801,000,000 0,005 0,010 0,015 0,020
0,025Constante de Torque% do maximo valor .%RPM%Potncia Grfico 2 -
Torque e RPM x Potncia
Paraobterosdadosdogrfico2etabela1,apenasotorquefoivariado.Para
constantesdetorquemenorque0,002emaiorque0,022omotorfuncionaporalgum
tempo e depois para. 44 Grfico 3 - Posio e Velocidade x Tempo
Simulando notempo de400 segundos, podemos observar no grfico 3 que
aps 100segundosavariaodevelocidadetemumvalormximoeumvalormnimo
estvel:-9rad/se-4rad/s,respectivamente.Ummelhorajustenocontrapeso,deve
diminuir essa variao. Essa estabilidade tambm percebida com relao
variao de temperatura do
fluidonosvolumesdecontrole.Nogrfico4,podemosobservaradiminuiona
temperaturadofluidonapartequente,umganhonapartefriaeumaconstnciana
parte intermediaria ao longo do
tempo.Conformepodeseranalisadopelasequaes3.12e3.13dametodologiade
termodinmicaetransfernciadecalor,avariaodetemperaturadofluidonoest
45 fortemente relacionada ao coeficiente de conveco. Quando maior
for o coeficiente de conveco menor ser a queda de temperatura.
Grfico 4 - Temperatura x Tempo
Damesmaformaquantomaioravelocidadederotaodomotor,maiorsera
velocidade das massas de ar que entram e saem, e menor ser o tempo
de contato do
fluidocomasuperfcie.Dessaformanohtemposuficienteparamantera
temperatura do fluido na mesma temperatura da superfcie. Esse
fenmeno pode ser observado no grfico 5, no qual foi feita uma
simulao com valor pequeno para a constante de torque (Kt), no caso
Kt=0,001. 46
Omotorganhamuitavelocidadenoinicio,poremastemperaturasdaspartesse
alteram drasticamente, fazendo o motor perder potncia e parar aps
algum tempo. Grfico 5 - Comportamento para Kt baixo
Notempoemtornode25segundosocorreuminversodosvaloresde temperatura.
Na parte da superfcie quente, a temperatura do fluido fica mais
fria que na
partedasuperfciefria.Mesmoassim,devidoinrcia,omotorcontinuafuncionando
at em torno de 40 segundos quando no consegue mais realizar uma
volta
completa.Atemperaturaaospoucosvaivoltandoaosvaloresiniciais,masnotemosmaisas
condies de posio e velocidade para o motor funcionar. 47 Voltando
para valores de Kt mais elevado (Kt=0,010), no qual ele fique em
regime permanente funcionando, conforme visto no grfico 3, vamos
agora observar no grfico 6 o comportamento da presso, fora e
temperatura. Grfico 6 - Comportamento da Presso, Fora e Temperatura
para regime permanente funcionando.
Vemosnogrficodapressoquetemospressesacimaeabaixodapresso
atmosfrica(1,01325x10-5Pa).Nocasoprticoapressopositiva,acimada
atmosfrica, tende a expulsar massa do motor. Na presso negativa,
abaixo da presso atmosfrica, o ar tende a entrar no motor.Podemos
notar que, tanto a Presso como a Fora, no tem simetria na oscilao.
Fica claro que a mxima fora positiva maior que a mxima fora
negativa. 48
Comrelaotemperaturamdiapodemosnotarquesuaamplitudediminuiat
atingiroregimepermanente,issoacontece,poisastemperaturasdaspartestambm
variam, conforme visto no grfico 4. Sendo a presso critica do ar
Pc=37,7*105 Pascal e a temperatura critica Tc=133K e que a mxima
presso p 1,7*105 Pascal e a mnima T 325K temos que a presso
reduzidamxima0,045eamnimatemperaturareduzida2,44.Dessaformade
acordo com a figura 3.12 de SHAPIRO[1] podemos tratar o ar desse
sistema no modelo ideal com resultados satisfatrios. Grfico 7 -
Velocidade e Posio em 25s
Nogrfico7podemosobservarqueat5segundosexisteumganhode velocidade,
aps isso ela vai diminuindo at se estabilizar em 100 segundos. 49
Grfico 8 - Condio inicial de 5rad/s positivo
Comopodemosvernogrficoacimaumavelocidadeinicialpositivade5rad/se
posio inicial de beta1 igual a zero radianos, o motor no funciona.
Foi simulado para
posioinicialigualaPiradianosevelocidadeinicialpositivaiguala5rad/seele
tambm no funcionou. 50 Grfico 9 - Condio Inicial de 15rad/s
positivo Aumentando a condio inicial de velocidade de beta1 para
15rad/s ele funciona, mas no sentido horrio. A posio inicial foi
com beta1 igual a zero.
Dessaformaanalisandoasgrficos3,7e9podemosafirmarqueomotortem rotao
nica de giro, no sentido horrio. 51 Grfico 10 - Volumes
Ovolumemorto,volumequenorealizatrabalho,muitoprejudicialparao
desempenho do motor [4]. Podemosobservar que aparte quente
apresenta quase 5ml (1 ml = 106 m) de volume morto. Como no estamos
alterando as dimenses, no foi possvel verificar o efeito do volume
morto. 52 Grfico 11 - Presso x Volume No ciclo ideal, em verde no
grfico 11, foi calculado apenas alguns pontos. A rea
internaotrabalhorealizado,dessaformapodemosnotaradiferenadocicloideal
para o ciclo simulado.O grfico 11 foi criado para o tempo variando
de 0 at 100 segundos. Como temos variao do trabalho realizado em
cada volta, o ciclo simulado tem vrios caminhos at atingir o regime
permanente, na rea interna menor. 53
ParaomotoroperandocomTqiguala600C(873K)fornece0,7556Watts,com um
toque de 0,0820N*m e RPM de 87,96. Nas mesmas condies mas com Tq
igual a 400C o motor fornece uma potencia de 0,653Watts. Dessa
forma um aumento em 50% na temperatura, temos um aumento de 15% na
potncia. 54 Concluso So inmeros os parmetros envolvidos no
desempenho do motor Stirling. Alguns
delesforamexemplificadosnessetrabalho.Ametodologiaempregadaapenasuma
tentativadesimulaodosefeitosqueocorremnaprticaecertamenteexistem
fenmenosquenoestoinclusonosresultados.Entretantovamosrealizaralgumas
anlises para somar na discusso do assunto. Como podemos ver no
grfico 1 uma pequena variao da quantidade do fluido de trabalho
altera significativamente o desempenho do motor. Na concepo do
projeto do
motormuitocomplicadoconstruirumacmeraquesemovimentaequenoaja
vazamento do fluido de trabalho, ainda mais com um torque to baixo
os ajustes entre o
pistoeocilindronopodesermuitojusto,queomotornovenceriaoatritoinicial.
Essa uma dificuldade tecnolgica.
interessanteobservarqueasensibilidadedomotorcomrelaoamassado
fluidodetrabalhotambmobservadocomrelaoaotorque.Comopodemosver
tantono grfico 2 e mais claramente no grfico 5, o torqueatingindo
certo valor, tanto
mximoquantomnimofazomotorperderpotnciadeformaqueelenotenha condies
de voltar ao regime permanente. Essa uma dificuldade de projeto. No
posso afirmar a viabilidade do motor Stirling por esse trabalho e
tambm por que motivo ele no largamente utilizado na sociedade, mas
talvez essas dificuldades
tecnolgicasedeprojetoresultemnaprticaemmotorespoucosfuncionaiseque
necessitemperiodicamenteajustarascondiesiniciaisparavoltarfuncionar.
Entretantonosltimosanosvemosinmerasempresassurgiremcomgrande
investimento produzindo o motor Sitling solar acoplado ao gerador
eltrico. interessante em trabalhos futuros analisar os efeitos da
variao das dimenses
daspeasbemcomooimpactodovolumemortoeposteriormenteconstruirum
prottipo para comparar os resultados. 55 Referncias Bibliogrficas
[1]MORAN,MichealJ.,SHAPIRO,HowardN.,PrincpiosdeTermodinmica para
Engenharia, John Wiley & Sons, 2003.
[2]INCROPERA,FrankP.,DEWITT,DavidP.,FundamentosdeTransferncia de
Calor e Massa, John Wiley & Sons, 2000.
[3]Listofmomentsofinertia.Wikipedia.Disponvelem:
http://en.wikipedia.org/wiki/List_of_moments_of_inertia
[4]DARLINGTON,Roy;STRONG,Keith.StirlingandHotAirEngines,2005.
Crowood Press. [5] DOUGHTY, S. Mechanics of Machines, John Wiley
& Sons, USA, 1988.
[6]BREGION,GregoryDaniel.AnliseDinmicadeumSistemaPino-Pisto com
Lubrificao Hidrodinmica. 2008. 56 ANEXO I Tabela das Propriedades
do AR como gs ideal Tabela 2 - Propriedade do AR retirado de [1] e
[2] Temperatura T [K] Cv Densidade [kg/m] Condutividade k[ W/(m*K)
] Entalpia H (KJ/Kg) 2000,7141,74580,0181199,97
2500,7161,39470,0223250,05 3000,7181,16140,0263300,19
3500,7210,99500,0300350,49 4000,7260,87110,0338400,98
4500,7330,77400,0373451,80 5000,7420,69640,0407503,02
5500,7530,63290,0439554,74 6000,7640,58040,0469607,02
6500,7760,53560,0497659,84 7000,7880,49750,0524713,27
7500,8000,46430,0549767,29 8000,8120,43540,0573821,95 57 ANEXO II
Clculo dos Coeficientes de Conveco
Novolumedecontroledoregeneradoroproblemadotipotuboconcntrico
anular. Conforme Incropera[2] temos coeficientes de conveco
distintos associados s superfcies interna e externa. Os nmeros de
Nusselt correspondentes so da forma: kD hNukD hNuh ssh ee== O
dimetro hidrulico Dh : e s hD D D =Para o caso de escoamento
laminar planamente desenvolvido com uma superfcie isolada e a outra
a uma temperatura constante os Nussels podem ser obtidos por tabela
(Tabela disponvel na pgina 346, Incropera[2]). Tanto no volume de
controle frio quanto no volume de controle quente idealizamos
acondiodetubocircularcomescoamentolaminarplenamentedesenvolvidocom
temperaturadesuperfcieuniforme.Nocasoestamossubestimandoonmerode
Nussels.Dessaformaaquedadetemperaturanosvolumesquenteeaaumentoda
temperaturanovolumefrioseromaisacentuadasnessasimulaoqueemumcaso
real. 58 ANEXO III Inrcia dos Elementos Tabela 3 - Momentos de
Inrcia, retirado deWikipedia[3] DescrioFiguraMomento de Inrcia
Cilindro Slido de raio r, altura h e massa m CuboSlidodealturah,
larguraw,comprimentode massa m. Pontodemassamcoma distncia r do
eixo de rotao I = mr2
Paraasbielasforamutilizadosainrciadecuboslido,assimforameliminados
todos os detalhes de cada pea. Como apenas o contra peso roda fora
do seu eixo de
centrodemassa,oteoremadoseixosparalelos,naterceiralinhadessatabelafoi
utilizado. Para o disco de inrcia foi utilizada a inrcia de
cilindro slido rodando em Z. 59 ANEXO IV Mecanismo de Sincronia:
calculo dos ngulos. Equaes de loop em x e em y: C M sen B sen B sen
BD M B B B= = + = = + 2 ) ( 1 ) ( 2 ) ( 34 ) cos( 1 ) cos( 2 ) cos(
31 2 31 2 3 Elevando ao quadrado temos: ( ) ( )( ) ( )22 3222322
32223) sin( ) sin( 2 3 2 ) ( 2 ) ( 3) cos( ) cos( 2 3 2 ) cos( 2 )
cos( 3C B B sen B sen BD B B B B= + + = + + Somando as equaes
temos: ( ) ( ) ( )() () ( )zz w awB BC D B BC D B B B BC D sen sen
B B sen B sen B == = = + + = + = + + + = + + + + + 2 33 22 2 2 23
22 23 22 22 22 3 2 3 2222 23232 2) cos(2 3 22 3) cos(cos( 2 3 2 1 2
1 3cos cos 2 3 2 cos 2 cos 3 Aplicando 3na equao de loop em x,
temos: ( )( )( ) ( ) ( )( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) )
cos( 2 cos 3 4 ) cos( 1 2 ) ( cos 2 cos 3 4 ) cos( 1 cos 1 ) ( 3)
cos( 2 cos 3 4 ) cos( 1 cos 1 ) ( 3) cos( 2 cos 3 4 ) cos( 1 cos 1
) ( 30 4 ) cos( 1 ) cos( 2 cos 1 ) ( 3 cos cos 30 4 ) cos( 1 ) cos(
2 ) ( cos cos 30 4 ) cos( 1 ) cos( 2 ) cos( 32 1 22 2 21 22 222 122
/ 1222 12 / 1221 22 / 12221 2 2 21 2 2 + + + = + =||
\| + = = + + + = + + + = + + B z B M B B z B M B z sen BB z B M
B z sen BB z B M B z sen BM B B z sen B z BM B B z sen sen z BM B B
z B Podemos arranjar a ultima equao comouma equaode segundo
graupara a varivel) cos(2 , assim temos: 60 ( ) ( ) [ ] ( ) ( ) [ ]
( ) [ ] ) ( 3 4 ) cos( 1 ) cos( 2 cos 3 4 ) cos( 1 2 ) ( cos 2 cos
3 ) ( 32 22 2 1 22 2 2z sen B M B B z B M B B z B z sen B + + + +
Dessa forma: c b a + + ) cos( ) ( cos2 22 Logo: ac a b b24) cos(22
= Portanto: |||
\| =ac a b ba24cos22 61 ANEXO V Programa em Matlab Arquivo
principal stirling.m %Stirling Tipo Beta - Analise do Mecanismo
clear all; close all; clc global t1 Tm Tq1 Vq1 Vq2 Tq2 guarda Vt1
Vt2 Vf1 Vf2 Vqf1 Vqf2 Tf1 Tf2 Tqf1 guarda2 AreaQ2 AreaQ1 AreaF2
AreaF1 Torque guarda3 Torque=0.01; %0.005 % [N*m/(rad/s)]
AreaQ2=AreaQ1; AreaF2=AreaF1; Tq1=400+273; Tf1=50+273; Tqf1=498;
Tm=464; % Temperatura Media inicial do Fluido Tq2=Tq1; Tf2=Tf1;
Tqf2=Tqf1; ti=0; tf=100; % n=500; % nmero de pontos de discretizao
% dt=(tf-ti)/n;% tamanho do passo da discretizao %% t=ti:dt:tf;
options = odeset('MaxStep',0.1); [T,Y]=ode45('lagrange',[ti
tf],[0;-5],options); % inicial =[0 ; -5] %Grafico figure
subplot(211) plot(T,Y(:,1)) title('Deslocamento - Beta1');
xlabel('tempo [s]'); ylabel('[rad]'); subplot(212) plot(T,Y(:,2))
title('Velocidade - Beta1'); xlabel('tempo [s]'); 62
ylabel('[rad/s]'); %%%%%%TESTE %----DADOS DE ENTRADA------------
B1=11; B2=34; B3=43; B4=78; B5=44; B6=12; B7=64; M1=24;
M2=(47-2.8086); M3=M2; %nao alterar, pois muda o mecanismo. M4=26;
M5=16; betaA=77*pi/180; Dc=35; %diametro interno da camisa/Diametro
do Pistao Cc=168; %comprimento da camisa Cd=95; %comprimento do
pistao deslocador Dd=33; % diametro pistao deslocador %Posiao dos
Centro de Massas uPb1=B1/2; vPb1=0; uPb4=B4/2; vPb4=0; %----LOOP DE
POSICAO----------- % beta1=Y(:,1); %Mecanismo de Potencia %em x
=> B1*cos(beta1) + B4*cos(beta4) - XP + Dp - M1 -M4 = 0 %em y
=> B1*sen(beta1) + B4*cos(beta4) = 0
beta4=asin(-B1*sin(beta1)/B4);
XPI=B1*cos(beta1)+B4*cos(beta4)+M5-M1-M4; % %Mecanismo de Sincronia
w=( (-M4+B1*cos(beta1)).^2 + (-M2+B1*sin(beta1)).^2 -B2.^2
-B3.^2)/(2*B3*B2); %equacao auxiliar z=acos(w); a=(B3*sin(z)).^2 +
(B3*cos(z)+B2).^2 ; b=-2*(B1*cos(beta1)-M4).*(B3*cos(z)+B2);
c=(B1*cos(beta1)-M4).^2 -(B3*sin(z)).^2;
beta2=-acos((-b+sqrt(b.^2-4.*a.*c))./(2*a)); beta3=beta2-acos(w);
63 beta5=betaA+beta3; % %Mecanismo de Deslocamento
beta6=asin(-M2-B5*sin(beta5))/B6;
XDE=real(B5*cos(beta5)+B6*cos(beta6)+B7-M1); %%%%
tamanho=size(Y,1);
RPM=((((Y(tamanho,1)-Y(1,1)))/(2*pi))/(tf-ti))*60
Tork=Torque*median(Y(:,2)) Potencia=Tork*2*pi*RPM/60
animar=input('Deseja animar: [1]-sim[4]-nao => '); if
(animar==1) xy=[0 0 B1*cos(beta1(1))B1*sin(beta1(1))
B1*cos(beta1(1))+B4*cos(beta4(1)) 0 M4 M2
B3*cos(beta3(1))+M4B3*sin(beta3(1))+M2 B5*cos(beta5(1))+M4
B5*sin(beta5(1))+M2 B5*cos(beta5(1))+M4+B6*cos(real(beta6(1)))0
B5*cos(beta5(1))+M4+B6*cos(real(beta6(1)))+B7 0]; a=[1 1 0 0 0 0 0
0 0 1 1 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 1 1 0 0 0 1 0 1 1 0 0 0 0
0 0 1 0 1 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 1 1]; figure
[xd,yd]=gplot(a,xy); cla; axis([-50 240 -90 190]); axis off;
set(gca,'Drawmode','Fast');
line([M4+M1 M4+M1+Cc M4+M1+Cc M4+M1],[Dc/2 Dc/2 -Dc/2 -Dc/2]);
patch([XPI(1)-M5+M4+M1 XPI(1)+M4+M1 XPI(1)+M4+M1 XPI(1)-M5+M4+M1
XPI(1)-M5+M4+M1],[Dc/2 Dc/2 -Dc/2 -Dc/2 Dc/2],'b'); 64
patch([XDE(1)+M4+M1 XDE(1)+M4+M1+Cd XDE(1)+M4+M1+Cd XDE(1)+M4+M1
XDE(1)+M4+M1],[Dd/2 Dd/2 -Dd/2 -Dd/2 Dd/2],'r');
M = moviein(tamanho);%te for i=1:tamanho, while i > tamanho
i=i-tamanho; end xy=[00 B1*cos(beta1(i)) B1*sin(beta1(i))
B1*cos(beta1(i))+B4*cos(beta4(i)) 0 M4 M2
B3*cos(beta3(i))+M4B3*sin(beta3(i))+M2
B5*cos(beta5(i))+M4B5*sin(beta5(i))+M2
B5*cos(beta5(i))+M4+B6*cos(real(beta6(i))) 0
B5*cos(beta5(i))+M4+B6*cos(real(beta6(i)))+B70];% gplot(a,xy);
axis([-50 240 -90 190]); axis off % % line([M4+M1 M4+M1+Cc M4+M1+Cc
M4+M1],[Dc/2 Dc/2 -Dc/2 -Dc/2]); patch([XPI(i)-M5+M4+M1
XPI(i)+M4+M1 XPI(i)+M4+M1 XPI(i)-M5+M4+M1 XPI(i)-M5+M4+M1],[Dc/2
Dc/2 -Dc/2 -Dc/2 Dc/2],'b'); patch([XDE(i)+M4+M1 XDE(i)+M4+M1+Cd
XDE(i)+M4+M1+Cd XDE(i)+M4+M1 XDE(i)+M4+M1],[Dd/2 Dd/2 -Dd/2 -Dd/2
Dd/2],'r'); %drawnow ; M(:,i) = getframe;end end % if animar %Para
animar novamente digite: movie(M) 65 Arquivo da equao diferencial
nomeado como lagrange.m function yp = lagrange(t,y) global t1 Tm
Tq1 Vq1 Vq2 Tq2 guarda Vt1 Vt2 Vf1 Vf2 Vqf1 Vqf2 Tf1 Tf2 Tqf1
guarda2 AreaQ2 AreaQ1 AreaF2 AreaF1 Torque guarda3 %Resolve a
equaao de movimento do Motor Stirling beta1 = y(1); %posicao
inicial do angulo beta1 %----DADOS DE ENTRADA------------ % %
%DIMENCOES DO MOTOR %DIMENCOES DAS BIELAS %Distancia Entre Eixos
B1=0.011; B2=0.034; B3=0.043; B4=0.078; B5=0.044; B6=0.012;
B7=0.064; %Largura lB1=0.003; lB2=0.003; lB3=0.003; lB4=0.003;
lB5=0.003; lB6=0.003; %Altura aB1=0.010; aB2=0.010; aB3=0.010;
aB4=0.010; aB5=0.010; aB6=0.010; %Comprimento; cB1=0.050;
cB2=0.047; cB3=0.055; cB4=0.095; cB5=0.055; cB6=0.025; %Posiao dos
Centro de Massas uPb1=B1/2; 66 vPb1=0; uPb2=B2/2; vPb2=0;
uPb3=B3/2; vPb3=0; uPb4=B4/2; vPb4=0; uPb5=-B5/2; vPb5=0;
uPb6=B6/2; vPb6=0; % %Dimensoes Extras %Disco de Inercia Ddi=0.092;
% Diametro Disco de Inercia %ContraPeso cCP=0.012; %Comprimento
aCP=0.007; %Altura lCP=0.003; %Largura uCP=-0.008; %Posicao centro
de massa em U vCP=0.013; %Posicao centro de massa em V %Dimenscoes
de Posiao e Tamanho M1=0.024; M2=(47-2.8086)*0.001; M3=M2; %nao
alterar, pois muda o mecanismo. M4=0.026; M5=0.016;
betaA=77*pi/180; Dc=0.035; %diametro interno da camisa/Diametro do
Pistao Cc=0.168; %comprimento da camisa Cd=0.095; %comprimento do
pistao deslocador Dd=0.033; % diametro pistao deslocador % % %DADOS
DO FLUIDO K_fluido=dado('k',Tm);% do AR [W/(m*K)] a
300Kcp_fluido=dado('cp',Tm); % do AR [kJ/(kg*K)] a 300K
RO_fluido=dado('ro',Tm); %kg/m^3 Pi=101325; %[Pascal] Pressao
Inicial do FluidoR=8.314472; %[J/(K*mol)]Constante Universal dos
Gases
n=0.000751;%Pi*((Cc*pi*Dc^2/4)-((Cd+0.072199)*pi*Dd^2/4))/(R*(298))
%%%%%%%ATENCAO colocar valor de T1 inica%%%% R*T1
massa_fluido=28.97*n; % [kg]de massa do fluido % % %DADOS DA
TEMPERATURA NO MOTOR Tsq=400+273; %Temperatura Parte Quente
[KELVIN] Tsf=50+273; %Temperatura Parte Fria [KELVIN]
Tsqf=(Tsq+Tsf)/2; %Temperatura da Parte do Meio [KELVIN] 67
%-------------% %RESOLVENDO % %-------------% %----LOOP DE
POSICAO----------- % %Mecanismo de Potencia %em x =>
B1*cos(beta1) + B4*cos(beta4) - XP + Dp - M1 -M4 = 0 %em y =>
B1*sen(beta1) + B4*cos(beta4) = 0 beta4=asin(-B1*sin(beta1)/B4);
XPI=B1*cos(beta1)+B4*cos(beta4)+M5-M1-M4; % %Mecanismo de Sincronia
w=( (-M4+B1*cos(beta1)).^2 + (-M2+B1*sin(beta1)).^2 -B2.^2
-B3.^2)/(2*B3*B2); %equacao auxiliar z=acos(w); a=(B3*sin(z)).^2 +
(B3*cos(z)+B2).^2 ; b=-2*(B1*cos(beta1)-M4).*(B3*cos(z)+B2);
c=(B1*cos(beta1)-M4).^2 -(B3*sin(z)).^2;
beta2=-acos((-b+sqrt(b.^2-4.*a.*c))./(2*a)); beta3=beta2-acos(w);
beta5=betaA+beta3; % %Mecanismo de Deslocamento beta6=asin(
(-M2-B5*sin(beta5))/B6); XDE=B5*cos(beta5)+B6*cos(beta6)+B7-M1;
%-----COEFICIENTES DE VELOCIDADE Ks=zeros(7,1);
Ks(1,1)=(B1/B2)*((cos(beta1)*tan(beta3)-sin(beta1))/(tan(beta3)*cos(beta2)-sin(beta2)));%Beta2
Ks(2,1)=(-B2/B3)*(cos(beta2)*Ks(1,1)+B1*cos(beta1))/cos(beta3);%Beta3
Ks(3,1)-(B1*cos(beta1)/(B4*cos(beta4)));%Beta4
Ks(4,1)=Ks(2,1);%Beta5
Ks(5,1)=(B5*cos(beta5)/(-B6*cos(beta6)))*Ks(4,1); %Beta6
Ks(6,1)=-B1*sin(beta1)-B4*sin(beta4)*Ks(3,1); %Xp
Ks(7,1)=-B5*sin(beta5)*Ks(4,1)-B6*sin(beta6)*Ks(5,1); %Xd
%----ACELERACAO L=zeros(7,1); % L(1,1)=(
-B3*cos(beta3)*Ks(2,1)^2-B2*cos(beta2)*Ks(1,1)^2+B1*cos(beta1)
)/(B2*sin(beta2)); %Beta2 L(2,1)=(
B3*sin(beta3)*Ks(2,1)^2+B2*sin(beta2)*Ks(1,1)^2-B1*sin(beta1)
)/(B3*sin(beta3));%Beta3
L(3,1)=(B1*sin(beta1)/(B4*cos(beta4)))+(B4*sin(beta4)/(B4*cos(beta4)))*Ks(3,1)^2;
%Beta4 L(4,1)= L(3,1); %Beta5 L(5,1)= (
B5*sin(beta5)*Ks(4,1)^2+B6*sin(beta6)*Ks(5,1)^2)/(B6*cos(beta6));
%Beta6 L(6,1)=-B1*cos(beta1)-B4*cos(beta4)*Ks(3,1)^2; %XPI L(7,1)=
-B5*cos(beta5)*Ks(4,1)^2-B6*cos(beta6)*Ks(5,1)^2; %Xd 68
%----COEFICIENTE DE VELOCIDADE---- Kc=zeros(24,1); %
Kc(1,1)=-uPb1*sin(beta1)-vPb1*cos(beta1); %XPB1
Kc(2,1)=uPb1*cos(beta1)-vPb1*sin(beta1); %YPB1
Kc(3,1)=(-uPb2*sin(beta2)-vPb2*cos(beta2))*Ks(1,1)-B3*sin(beta3)*Ks(2,1);
%XPB2
Kc(4,1)=(uPb2*sin(beta2)-vPb2*sin(beta2))*Ks(1,1)+B3*cos(beta3)*Ks(2,1);
%YPB2 Kc(5,1)=(uPb3*sin(beta3)+vPb3*cos(beta3))*Ks(2,1); %XPB3
Kc(6,1)=(-uPb3*cos(beta3)+vPb3*sin(beta3))*Ks(2,1); %YPB3
Kc(7,1)=-B1*sin(beta1)-(uPb4*sin(beta4)+vPb4*cos(beta4))*Ks(3,1);
%XPB4
Kc(8,1)=B1*cos(beta1)+(uPb4*cos(beta4)-vPb4*sin(beta4))*Ks(3,1);
%YPB4 Kc(9,1)=-uPb5*sin(beta5)*Ks(4,1)-vPb5*cos(beta5)*Ks(4,1);
%XPB5 Kc(10,1)=uPb5*cos(beta5)*Ks(4,1)-vPb5*sin(beta5)*Ks(4,1);
%YPB5
Kc(11,1)=-B5*sin(beta5)*Ks(4,1)-uPb6*sin(beta6)*Ks(5,1)-vPb6*cos(beta6)*Ks(5,1);
%XPB6
Kc(12,1)=B5*cos(beta5)*Ks(4,1)+uPb6*cos(beta6)*Ks(5,1)-vPb6*sin(beta6)*Ks(5,1);
%YPB6 Kc(13,1)=-B1*sin(beta1)-B4*sin(beta4)*Ks(3,1); %XPP
Kc(14,1)=-B5*sin(beta5)*Ks(4,1)+B6*cos(beta6)*Ks(5,1); %XPD
Kc(15,1)=1; %Beta1 Kc(16,1)=Ks(1,1); %Beta2 Kc(17,1)=Ks(2,1);
%Beta3 Kc(18,1)=Ks(3,1); %Beta4 Kc(19,1)=Ks(4,1); %Beta5
Kc(20,1)=Ks(5,1); %Beta6 Kc(21,1)=-uCP*sin(beta1)-vCP*cos(beta1);
%xCP Kc(22,1)=uCP*cos(beta1)-vCP*sin(beta1);%yCP Kc(23,1)=1;
%Inercia CP Kc(24,1)=1; %Disco de Inercia %----COEFICIENTE DE
ACELERACAO---- LA=zeros(24,1); %
LA(1,1)=-uPb1*cos(beta1)+vPb1*sin(beta1); %XPB1
LA(2,1)=-uPb1*sin(beta1)-vPb1*cos(beta1); %%YPB1
LA(3,1)=(-uPb2*cos(beta2)+vPb2*sin(beta2))*Ks(1,1)^2-B3*cos(beta3)*Ks(2,1)^2;
%XBP2 LA(4,1)=(
uPb2*cos(beta2)-vPb2*cos(beta2))*Ks(1,1)^2-B3*sin(beta3)*Ks(2,1)^2;
%YPB2 LA(5,1)=(uPb3*cos(beta3)-vPb3*sin(beta3))*Ks(2,1)^2; %XPB3
LA(6,1)=(uPb3*sin(beta3)+vPb3*cos(beta3))*Ks(2,1)^2; %YPB3
LA(7,1)=-B1*cos(beta1)-(uPb4*cos(beta4)+vPb4*cos(beta4))*Ks(3,1)^2;
%XPB4
LA(8,1)=-B1*sin(beta1)-(uPb4*sin(beta4)+vPb4*cos(beta4))*Ks(3,1)^2;
%YPB4 LA(9,1)=-uPb5*cos(beta5)*Ks(4,1)^2+vPb5*sin(beta5)*Ks(4,1)^2;
%XPB5
LA(10,1)=-uPb5*sin(beta5)*Ks(4,1)^2-vPb5*cos(beta5)*Ks(4,1)^2;
%YPB5
LA(11,1)=-B5*cos(beta5)*Ks(4,1)^2-uPb6*cos(beta6)*Ks(5,1)^2-vPb6*sin(beta6)*Ks(5,1)^2;
%XPB6
LA(12,1)=-B5*sin(beta5)*Ks(4,1)^2-uPb6*sin(beta6)*Ks(5,1)^2-vPb6*cos(beta6)*Ks(4,1)^2;
%YPB6 LA(13,1)=-B1*cos(beta1)-B4*cos(beta4)*Ks(3,1)^2; %XPP
LA(14,1)=-B5*cos(beta5)*Ks(4,1)^2-B6*sin(beta6)*Ks(5,1)^2; %XPD
LA(15,1)=0; %Beta1 LA(16,1)=L(1,1); %Beta2 LA(17,1)=L(2,1); %Beta3
LA(18,1)=L(3,1); %Beta4 LA(19,1)=L(4,1); %Beta5 LA(20,1)=L(5,1);
%Beta6 LA(21,1)=-uCP*cos(beta1)+vCP*sin(beta1); %xCP 69
LA(22,1)=-vCP*sin(beta1)-vCP*cos(beta1); %yCP LA(23,1)=0; %Inercia
CP LA(24,1)=0; %Disco de Inercia %----MATRIZ DE INERCIA--------
M=zeros(24,24); % M(1,1)=0.3; %massa B1 M(3,3)=0.3; %massa B2
M(5,5)=0.3; %massa B3 M(7,7)=0.3; %massa B4 M(9,9)=0.3; %massa B5
M(11,11)=0.3; %massa B6 M(13,13)=0.4; %massa Pp M(14,14)=0.4;
%massa Pd M(21,21)=0.1765; %massa ContraPeso(CP) massaDI=3; %massa
do disco de inercia M(2,2)=M(1,1); M(4,4)=M(3,3); M(6,6)=M(5,5);
M(8,8)=M(7,7); M(10,10)=M(9,9); M(12,12)=M(11,11);
M(22,22)=M(21,21); M(15,15)=(1/12)*M(1,1)*(cB1^2+aB1^2); %Inercia
B1 (1/12)*M(1,1)*(B1^2+aB1^2)+M(1,1)*(uPb1)^2;
M(16,16)=(1/12)*M(3,3)*(cB2^2+aB2^2); %Inercia B2
M(17,17)=(1/12)*M(5,5)*(cB3^2+aB3^2)+M(5,5)*(uPb3)^2; %Inercia B3
M(18,18)=(M(7,7)/12)*(cB4^2+aB4^2); %Inercia
B4M(19,19)=(1/12)*M(9,9)*(cB5^2+aB5^2)+M(9,9)*(cB5/2)^2; %Inercia
B5 M(20,20)=(1/12)*M(11,11)*(cB6^2+aB6^2); %Inercia B6
M(23,23)=(1/12)*M(21,21)*(cCP^2+aCP^2)+M(21,21)*(uCP^2+vCP^2);
%Inercia CP M(24,24)=massaDI*Ddi^2/2; %Inercia Disco de Inercia
%----MATRIZ AUXILIAR-------------- N1=0.5*(LA'*M*Kc+Kc'*M*LA);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %----VOLUMES
DE TRABALHO------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %VOLUMES
Vf1=(XDE-XPI)*pi*((Dc)^2)/4; %m^3 Vq1=(Cc-(XDE+Cd))*pi*((Dc)^2)/4;
%m^3 Vqf1=pi*(((Dc-Dd))^2)/4*Cd; %m^3 Vt1=Vf1+Vqf1+Vq1; %m^3 70
%AREA AreaF1=(XDE-XPI)*pi*((Dc))+(Dc^2)*pi/4;
AreaQ1=(Cc-(XDE+Cd))*pi*((Dc))+(Dc^2)*pi/4; AreaQFe1=Cd*pi*Dc;
AreaQFi1=Cd*pi*Dd; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Passo da resolucao da
funcao ODE %%% t1=[t1 t]; passo=size(t1,2); if passo == 1
InterTEMPO=t; Vf2=Vf1; Vqf2=Vqf1; Vq2=Vq1; Vt2=Vt1; Tq2=Tq1;
Tf2=Tf1; Tqf2=Tqf1; else InterTEMPO=t-t1(passo-1); end %% %Pressao
% %---CALCULO DA PROXIMA TEMPERTURA MEDIA------ %Calculo do
Coediciente de Conveccao NuFQ=3.66; % laminar planamente
desenvolvido, Ts uniforme, Pr>=0.6 [Pg.350,INCROPERA]
h_Q=NuFQ*dado('k',Tq1)/(Dc); h_F=NuFQ*dado('k',Tf1)/(Dc);
Q_F=h_F*(Tsf-Tf1)*AreaF1; %Watts Q_Q=h_Q*(Tsq-Tq1)*AreaQ1; %Watts
[NuMe,NuMs]=NuTubAnu(Dd/Dc); % [Nue,Nus] => Nu de Tubo
Concentrico Anular Dh=Dc-Dd; h_Me=NuMe*dado('k',Tm)/(Dh);
h_Ms=NuMs*dado('k',Tm)/(Dh); Q_Me=h_Me*(Tsqf-Tqf1)*AreaQFe1; %Watts
71 Q_Ms=h_Ms*(Tsqf-Tqf1)*AreaQFi1; %Watts Q_M=Q_Me+Q_Ms;
%TEMPERATURA MEDIA
Tm=(Tq1*Tf1*Tqf1)*Vt1/(Tf1*Tqf1*Vq1+Tq1*Tqf1*Vf1+Tq1*Tf1*Vqf1);
P=n*R*Tm/Vt1; %Pressao Forca=(P-Pi)*pi*(Dc)^2/4; %massas sendo
transferiada if (Vq2-Vq1) < 0m1=dado('ro',Tq1)*(Vq2-Vq1);%massa
saiu da parte quente Wq=0; Em1=m1*dado('h',Tq1);
m2=dado('ro',Tqf1)*(Vf2-Vf1);%massa que entra na parte fria
Wf=-P*(Vt1-Vt2); Em2=m2*dado('h',Tqf1);
elsem1=dado('ro',Tqf1)*(Vq2-Vq1); %massa que entra na parte quente
Wq=P*(Vt2-Vt1); Em1=m1*dado('h',Tqf1); m2=dado('ro',Tf1)*(Vf2-Vf1);
%massa que sai na parte fria Wf=0; Em2=m2*dado('h',Tf1); end
%BALANCO DE ENERGIA % %no volume de controle quente Tq2=(
(Q_Q*InterTEMPO-Wq+Em1)/dado('cv',Tq1) ) + Tq1; % %no volume de
controle frio Tf2=( (Q_F*InterTEMPO-Wf+Em2)/dado('cv',Tf1) ) + Tf1;
% %no volume de controle medio Tqf2=(Q_M*InterTEMPO+(-Em2-Em1)
)/dado('cv',Tqf1)+Tqf1; %atualizando as variaveis para o passo
seguinte Tq1=Tq2; Tf1=Tf2; Tqf1=Tqf2; 72 Vq2=Vq1; Vf2=Vf1; Vt2=Vt1;
Vqf2=Vqf1; AreaQ2=AreaQ1; AreaF2=AreaF1; guarda=[guarda Vt1];
guarda2=[guarda2 P]; %guarda3=[guarda3 Tqf1];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %----FORMULACAO DE LAGRANGE----
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %-----Equaes de primeira ordem
--------- % % ---------------------------FORMULACAO DE
LAGRANGE----------------------------- % %{ } {} { }{} % [Mg] *
{qpp} + 2 * (qp*[N1]) * {qp} - {{qp}'*[N1]*{qp}}= {Qf} %{ } {} {
}{} Mg = Kc'*M*Kc; % matriz de massa generalizada
Qf=Forca*Kc(13,1); %M1%M2 %M3%M4 %M5 %M6%CP dEPdq=(
M(1,1)*Kc(2,1)*Kc(15,1)+M(3,3)*Kc(4,1)*Kc(16,1)+M(5,5)*Kc(6,1)*Kc(17,1)+M(7,7)*Kc(8,1)*Kc(18,1)+M(9,9)*Kc(10,1)*Kc(19,1)+M(11,11)*Kc(12,1)*Kc(20,1)+M(21,21)*Kc(22,1)*Kc(15,1)
)*9.81; %Variacao da Energia Potencial %%% %%%EQUACOES %%%
yp(1,1)=y(2); yp(2,1)=inv(Mg)*(-y(2)*N1*y(2)-dEPdq
-Torque*y(2)+Qf); 73
Arquivosauxiliaresdefunodeapoio.NuTubAnu.mparacalcularonmerode
Nussel no tubo anular. function [Nue,Nus] = NuTubAnu(taxa) %Retorna
o Nu (numero de Nussel) para um tubo anular % % [Nue,Nus] =
NuTubAnu(taxa) % %Nue=Nussel Interno %Nus=NUssel da Superficie
%taxa= tNUe=[20.0017.4611.567.375.744.86]; tNUs=[ 3.66 4.06
4.114.234.434.86]; tD=[ 00.050.10.25 0.5 1];
Nue=interp1(tD,tNUe,taxa); Nus=interp1(tD,tNUs,taxa); Para
encontrar as propriedades do ar com o nome de dado.m function valor
= dado(v_retorno,valor,v_busca); %PROPIEDADES DO AR%
%dado(retorna,valor,busca) % %EXEMPLOS: %%dado('u',300) =
214.0700Retornou o valor de 'u' para o valor de 'te' igual a 300K %
%dado('te',214.07,'u') = 300 Retornou o valor de 'te' para o valor
de 'u' igual a 214.07 % % parametro => 'u'= energia interna
[kJ/kg] % parametro => 'ro' = densidade [kg/m^3] % parametro
=> 'cp' = calor especifico [kJ/(kg*K)] % parametro => 'k'=
condutividade [W/(m*K)] % parametro => 'te' = temperatura [K] %
parametro => 'h'= entalpia [kJ/kg] temp=[ 200250
300350400450500550600650700750800 850 900 950 1000]; % temperatura
[K] u=[ 142.56 178.28 214.07 250.02 286.16 332.62 359.49 396.86
434.78 473.25 512.33 551.99 592.30 633.18 674.58 716.55 758.94]; %
energia interna [kJ/kg] 74 ro=[1.7458 1.3947 1.1614 0.9950 0.8711
0.7740 0.6964 0.6329 0.5804 0.5356 0.4975 0.4643 0.4354 0.4097
0.3868 0.3666 0.3482]; %densidade [kg/m^3]
cp=[1.0071.0061.0071.0091.0141.0211.0301.0401.0511.0631.0751.0871.099
1.1101.1211.1311.141];k=[0.0181 0.0223 0.02630.030 0.0338 0.0373
0.0407 0.0439 0.0469 0.0497 0.0524 0.0549 0.05730.05960.0620.0643
0.0667]; % condutividade [W/(m*K)] h=[199.97 250.05 300.19 350.49
400.98 451.80 503.02 554.74 607.02 659.84 713.27 767.29
821.95877.18 932.93 989.24 1046.04]; cv=[0.714
0.7160.7180.7210.7260.7330.7420.7530.7640.7760.7880.8000.812
0.8230.8340.8445 0.855]; if v_retorno=='u' parametro=u; elseif
v_retorno=='ro' parametro=ro; elseif v_retorno=='cp' parametro=cp;
elseif v_retorno=='k' parametro=k; elseif v_retorno=='te'
parametro=temp; elseif v_retorno=='h' parametro=h; elseif
v_retorno=='cv' parametro=cv; end %verifica se foi dado o terceiro
parametro if nargin < 3 v_busca='te'; end if v_busca=='te'
busca=temp; elseif v_busca=='u' busca=u; elseif v_busca=='ro'
busca=ro; end
valor=interp1(busca,parametro,valor);
