UNESP FACULDADE DE ENGENHARIA DO CAMPUS DE GUARATINGUETÁ GUARATINGUETÁ 2012 2 GENEROSO NIEDERAUER DE OLIVEIRA MODELO DINÂMICO DE UM TÚNEL DE VENTO DE SOPRO Tese apresentada à Faculdade de Engenharia do Campus de Guaratinguetá, Universidade Estadual Paulista, para a obtenção do título de Doutor em Engenharia Mecânica na área de Transmissão e Conversão de Energia. Orientador: Prof. Dr. Luiz Roberto Carrocci Co-Orientador: Prof. Dr. Maurício Guimarães da Silva Co-Orientador: Prof. Dr. Maurício Araújo Zanardi Guaratinguetá 2012 3 Oliveira, Generoso Niederauer O48t Modelo Dinâmico de um Túnel de Vento de Sopro/ Generoso Niederauer de Olivera. - Guaratinguetá : [s.n.], 2012 136 f.: il. Referências: f. 111-114 Tese (Doutorado) – Universidade Estadual Paulista, Faculdade de Engenharia de Guaratinguetá, 2012 Orientador: Prof. Dr. Luiz Roberto Carrocci Co orientador: Prof. Dr. Maurício Guimarães da Silva Co orientador: Prof. Dr. Maurício Araújo Zanardi 1. Mecânica do fluídos I. Título CDU 532 4 UNESP UNIVERSIDADE ESTADUAL PAULISTA Faculdade de Engenharia do Campus de Guaratinguetá "MODELO DINÂMICO DE UM TÚNEL DE VENTO DE SOPRO" GENEROSO NIEDERAUER DE OLIVEIRA ESTA TESE FOI JULGADA ADEQUADA PARA A OBTENÇÃO DO TÍTULO DE “DOUTOR EM ENGENHARIA MECÂNICA” PROGRAMA: ENGENHARIA MECÂNICA ÁREA: TRANSMISSÃO E CONVERSÃO DE ENERGIA APROVADA EM SUA FORMA FINAL PELO PROGRAMA DE PÓS-GRADUAÇÃO Prof. Dr. José Antônio Perrella Balestieri Júnior Coordenador BANCA EXAMINADORA: Prof. Dr. LUIZ ROBERTO CARROCCI – DEN/FEG Prof. Dr. VICTOR ORLANDO GAMARRA ROSADO – DME/FEG Prof. Dra. VALESCA ALVES CORRÊA – UNITAU Prof. Dr. PEDRO JOSÉ DE OLIVEIRA NETO – DCTA/IAE/ALA Prof. Dr. LUIS FERNANDO GOUVEIA DE MORAES – DCTA/IAE/ALA Agosto de 2012 5 DADOS CURRICULARES GENEROSO NIEDERAUER DE OLIVEIRA NASCIMENTO 14.12.1943 – SANTA MARIA / RS FILIAÇÃO Generoso Zaballa de Oliveira Nilza Niederauer de Oliveira 1963/1967 Curso de Graduação Engenharia Mecânica - Universidade Federal de Rio Grande – RS. 2002/2004 Curso de Especialização em Certificação de Aeronaves Civis – Instituto Tecnológico de Aeronáutica – ITA e IFI – DCTA –Depto de Ciência e Tecnologia Aeroespacial. 2002/2005 Mestrado Profissionalizante - Instituto Tecnológico de Aeronáutica – ITA. 6 AGRADECIMENTOS Agradeço a Deus, a Jesus e aos queridos amigos de outras dimensões da vida. Agradeço muito aos familiares e demais entes queridos que me incentivaram incansavelmente. Também agradeço muito aos Prof. Dr. Luiz Roberto Carrocci e Prof. Dr. Maurício Araújo Zanardi que me cativaram profundamente como seres humanos que buscam de forma cristã e exemplar desempenhar suas missões de professores jamais desanimando ante qualquer dificuldade. Agradeço muito a atenção e amparo do Prof. Dr. Maurício Guimarães da Silva, ser humano de grande sensibilidade e percepção que tem a rara capacidade de transformar aparentes fracassos em vitórias inesquecíveis. Agradeço muito também a Sra. Maria Luiza da Pós-Graduação FEG, na área de Conversão e Transmissão de Energia, que com muita paciência e competência sempre conseguiu encaminhar/confirmar o recebimento com sucesso dos meus e- mails aos respectivos destinatários. 7 OLIVEIRA, G. N. Modelo Dinâmico de um Túnel de Vento de Sopro. 2012. xxxf. Tese (Doutorado em Engenharia Mecânica) – Faculdade de Engenharia do Campus de Guaratinguetá, Universidade Estadual Paulista, Guaratinguetá, 2012. RESUMO Neste trabalho é desenvolvido e apresentado anteprojeto de Túnel de Vento de Sopro (TVS). Para esta finalidade foi desenvolvida ferramenta computacional para análise aerotermodinâmica para Túneis de Vento de Sopro. Neste âmbito é apresentado o estado da arte dos TVS, respectiva descrição sumária de TVS, principais problemas operacionais, pré-dimensionamento com base nos requisitos operacionais, formulação matemática do túnel e do respectivo sistema de controle. Também é feita a apresentação das ferramentas computacionais com base nos requisitos que possam validar a metodologia proposta e são apresentados resultados e respectivas análises destes resultados. Na parte final deste trabalho é apresentada a respectiva validação e os principais resultados obtidos com um anteprojeto de túnel acadêmico. PALAVRAS-CHAVE: Mecânica Dos Fluidos / Túneis de vento de sopro / Projeto / Escoamento Supersônico / Ondas De Choque/ Escoamento Transiente Compressível 8 OLIVEIRA, G.N. Dynamic Model for a Blowdown Wind Tunnel. 2012. xxx pages. Thesis (Doctorate in Mechanical Engineering) - Faculdade de Engenharia do Campus de Guaratinguetá, Universidade Estadual Paulista, Guaratinguetá, 2012. ABSTRACT In this work it is developed a basic design of a Supersonic Blowdown Wind Tunnel (SBDWT). For this purpose computer programs were developed for design and thermo-aerodynamics analysis for a SBDWT. In this framework it is presented the state of art of SBDWT, brief description of a SBDWT, its main operational problems, pre-design based on operational requirements, and tunnel’s mathematical formulation. Also it is described the computational programs based on tunnel’s requirements; results are presented, and the corresponding analysis is also discussed. At the last part of this work it is presented the assumptions made for the mathematical modeling, the main results obtained and the validation of the computer programs are emphasized, as well. A the last part of this work there is presented the respective validation and the main results got for an academic draft design of a wind tunnel. KEYWORDS: Fluid Mechanics / Blowdown Supersonic Wind Tunnel / Design / Supersonic Flow / Shock Waves / Compressible Transient Flow 9 LISTA DE FIGURAS Figura 1. 1- Desenho esquemático de um túnel de vento de sopro ......................................... 22 Figura 2. 1 - Desenho esquemático de um túnel de vento de Sopro ....................................... 28 Figura 2. 2 - Fotografia Schlirien de um difusor supersônico (Neumann, 1949) .................. 34 Figura 2. 3 – Exemplos de Difusor Supersônico (Ref.: Neumann, Lustwerk, 1949) ............ 35 Figura 2. 4 – Espaço de Estados - Diagrama de Blocos .......................................................... 43 Figura 2. 5 - Diagrama de blocos em malha fechada para o TVS .......................................... 46 Figura 2. 6 - Diagrama SIMULINK: Controlador de Pressão de Estagnação ..................... 51 Figura 3. 1 – Geometria do problema, MOECKEL (1949) .................................................... 57 Figura 4. 1 - Razão de compressão em túneis de vento de Sopro ........................................... 71 Figura 4. 2 – Componente Longitudinal da Turbulência. ...................................................... 75 Figura 4. 3 – Componente Lateral da Turbulência. ............................................................... 76 Figura 4. 4 – Modelos de difusor com geometria fixa ............................................................. 79 Figura 4. 5 – Análise de difusores ............................................................................................ 80 Figura 5. 1– Simulação sem controlador (Mach= 3) ............................................................... 85 Figura 5. 2– Controle de 0P – Mach = 3,0 .............................................................................. 87 Figura 5. 3– Controle de STRe – Mach = 3,0 ......................................................................... 89 Figura 5. 4– Análise comparativa: TT – Mach = 3,0 .............................................................. 91 Figura 5. 5– Análise comparativa: abertura da válvula controle – Mach = 3,0 .................. 92 Figura 5. 6 – Análise comparativa: Desempenho do difusor – Mach = 3,0 ; 98,0�� ........ 94 Figura 5. 7 – Análise comparativa: Desempenho do difusor – Mach = 3,0 ; � = 0,68 ..... 95 Figura 5. 8 – Túnel de vento projetado: controle de 0P .......................................................... 96 Figura 5. 9 – Túnel de vento projetado: variação do ângulo de abertura ............................. 97 Figura 5. 10 – Túnel de vento projetado: número de Reynolds na seção de testes .............. 97 Figura 5. 11 – Túnel de vento projetado: análise de controladores ....................................... 100 Figura 5. 12 – Geometria do modelo em estudo....................................................................... 101 Figura 5. 13 – Distribuição de número de Mach na parede do modelo ................................. 102 Figura 5. 14 – Análise comparativa: forma e posição da onda de choque (Mach= 1,5) ...... 103 Figura 5. 15 – Erro relativo verificado no algoritmo de busca (Mach= 1,5) ......................... 103 Figura 5. 16 – Forma da onda de choque (Mach=1,5 ; 15� =15%) ........................................ 104 10 Figura 5. 17 – Identificação de k (Mach = 1,5 ; 15� =15%) .................................................... 105 Figura 5. 18 – Forma da onda de choque (Mach=3,0 ; 10� =10%) ........................................ 106 Figura 5. 19 – Identificação de k (Mach = 3,0 ; 10� =10%) .................................................... 106 11 LISTA DE TABELAS Tabela 2. 1– Relação entre gC e � ........................................................................................... 31 Tabela 2. 2– Equações linearizadas: Reservatório de ar e Câmara de tranquilização ........ 42 Tabela 2. 3– Modelo matemático para o TVS .......................................................................... 50 Tabela 2. 4– Entrada de Dados ................................................................................................. 51 Tabela 3. 1– Versões do código Missile Datcom ...................................................................... 55 Tabela 4. 1 – Requisitos para concepção do modelo .............................................................. 66 Tabela 4. 2 – Requisitos para concepção do TVS ................................................................... 67 Tabela 4. 3 – Exemplo de dimensionamento de um TVS ........................................................ 81 Tabela 5. 1– Organização da apresentação de resultados: Túnel de Vento ......................... 82 Tabela 5. 2– Organização da apresentação de resultados: Seção de Testes ......................... 83 Tabela 5. 3– Condições iniciais: FUNG(1987)......................................................................... 84 Tabela 5. 4– Quadro comparativo entre resultados experimentais e numéricos................. 87 Tabela 5. 5 – Condições de simulação ...................................................................................... 102 12 LISTA DE SÍMBOLOS m� Vazão em Massa Kg/s A Área m2 CD Coeficiente de descarga do bocal - Cg “Coeficiente de Projeto para Gás” da Válvula de - cp Calor Específico do Ar a Pressão Constante J/kg-K cv Calor Específico do Ar a Volume Constante J/kg-K D Diâmetro/Comprimento m d Menor Dimensão do Modelo m E(s) Erro do sinal no domínio da transformada de Laplace - e(t) Erro do sinal no domínio do tempo - G Vazão mássica por unidade de área (kg/s)/m2 G(s) Função de Transferência do Controlador PI - h Entalpia Específica J/kg H Altura do modelo m Ki Ganho Integral do Controlador - Kp Ganho Proporcional do Controlador - L Maior Dimensão do Modelo m l Dimensão/comprimento de referência do modelo. m M Número de Mach - n Coeficiente politrópico - P Pressão Total Pa p Pressão Estática Pa Q Vazão em volume m3/s q Transferência de Calor W/m2 R Raio (m) ou Const. dos Gases Perfeitos para o Ar = 287 J/kg-K Re, Rey Número de Reynolds - s Variável Transformada de Laplace - T Temperatura K t Tempo s 13 U Energia J V Volume m3 V0 Volume da câmara de Tranquilização m3 v Velocidade do Ar m/s a Velocidade do som = TR ��� m/s W Largura do modelo m B Componente da eq. da viscosidade de Sutherland Ns/m2K0, ou identificação de tela - C Componente da eq. da viscosidade de Sutherland K K0 Coef. de perda de pressão na tela - KL Coef. de perda de pressão em cada tela - fu Fator de redução do nível de turbul. axial - fv Fator de redução do nível de turbul. transversal - LISTA DE SÍMBOLOS - CARACTERES GREGOS � Porosidade da tela m2 Espaçamento entre os fios da tela m � Solidez da tela m2 Constante na equação (1) - � Diferença entre dois valores da mesma entidade física - � Massa específica do Ar kg/m3 �(t) Ângulo do Obturador da Válv. de Control., no tempo rad./graus (s) Ângulo do Obturador da Válvula Controladora no campo das transformadas de Laplace. - � Relação entre os Calores Específicos do Ar - � Viscosidade do ar kg�m/s � Máxima perda de pressão na válv. controladora - � Eficiência % � Razão de pressões = pressão considerada/pressão - 14 LISTA DE SÍMBOLOS – SUBSCRITOS * Condição Sônica ou garganta do Bocal de Aceleração - 0 Estagnação, Câmara de Tranquilização - 0s Saída do difusor - a Ar - Amb Ambiente - atm Atmosférica - BA Bocal de Aceleração D Difusor, Descarga - d Difusor - difusor Difusor - e Entrada - estag Estagnação - G Garganta - G BA Garganta do Bocal de Aceleração - G Garganta do Difusor - Inicial Inicial - m Médio, modelo - max Máximo ou Máxima - min Mínimo ou Mínima - rad Radianos - run Corrida s Choque - s, ed Saída do Difusor - sd Posição do choque no difusor - ST, se Seção de Testes, seção de ensaios - T Tanque de Ar - t Garganta (“throat”) do bocal de aceler. - td Garganta do Difusor - v Válvula de Controle - x, 0x, Pos. do choque, imediatam. antes do choque - y, 0y, Pos. do choque, imediatam. após o choque - 1 Estado 1 - 2 Estado 2 - 15 CTR Câmara de Tranquilização - LISTA DE SÍMBOLOS SUPERSCRITOS f Final - i Inicial - Set Valor objetivado - x Posição do choque ou imediatam. antes do choque - y Posição do choque ou imediatam. após o choque - n Coeficiente politrópico ou qtd. de telas em série - 16 SUMÁRIO 1. INTRODUÇÃO ....................................................................................................... 18 1.1 Definição do Problema ................................................................................................... 18 1.2 Objetivo .......................................................................................................................... 20 1.3 Túnel de Vento de Sopro ................................................................................................ 20 1.4 Revisão bibliográfica ...................................................................................................... 22 1.5 Escopo do Trabalho ........................................................................................................ 26 2. FORMULAÇÃO MATEMÁTICA PARA O TVS .............................................. 28 2.1 Hipóteses ........................................................................................................................ 28 2.2 Reservatório de Ar Comprimido .................................................................................... 29 2.3 Câmara de Tranquilização .............................................................................................. 31 2.4 Bocal do TVS ................................................................................................................. 32 2.5 Seção de Testes e Difusor .............................................................................................. 33 2.6 Sistema de Controle........................................................................................................ 37 2.7 Controlador ..................................................................................................................... 38 2.7.1 Controlador baseado na pressão de estagnação da câmara de Tranquilização ................. 39 2.7.2 Controlador baseado no número de Reynolds da seção de testes ..................................... 39 2.7.3 Determinação dos ganhos do controlador ......................................................................... 41 2.8 Implementação Numérica ............................................................................................... 47 2.9 Metodologia de integração ............................................................................................. 48 3. SEÇÃO DE TESTES .............................................................................................. 52 3.1 Missile Datcom............................................................................................................... 53 3.2 Onda de choque destacada ............................................................................................. 54 3.3 Formulação Matemática ................................................................................................. 57 3.4 Geometria da onda de choque ........................................................................................ 58 3.5 Distância entre o choque e o nariz do corpo .................................................................. 60 3.6 Princípio da continuidade ............................................................................................... 61 3.7 Algoritmo Numérico ...................................................................................................... 62 3.8 Método de Identificação de Sistemas ............................................................................. 64 4. DIMENSIONAMENTO DE UM TVS .................................................................. 66 4.1 Requisitos ....................................................................................................................... 66 4.2 Bocal de Aceleração ....................................................................................................... 68 4.3 Reservatório de Ar Comprimido .................................................................................... 68 4.4 Câmara de Tranquilização .............................................................................................. 72 4.5 Seção de Testes .............................................................................................................. 76 4.6 Difusor ............................................................................................................................ 77 4.7 Dimensionamento preliminar de um TVS...................................................................... 81 5 RESULTADOS ....................................................................................................... 82 5.1 Introdução ....................................................................................................................... 82 5.2 Validação da Metodologia de Anteprojeto de TVS ....................................................... 83 5.2.1 Condições Iniciais .................................................................................................... 83 5.2.2 Simulação de TVS Sem Controlador .............................................................................. 84 5.2.3imulação de TVS com Controle: Pressão de Estagnação .................................................. 85 5.2.4 Simulação de TVS com Controle: Número de Reynolds .............................................. 88 17 5.3 Análise Comparativa entre as Metodologias de Controle .............................................. 89 5.3.1 Temperatura de estagnação no reservatório de ar comprimido ..................................... 90 5.3.2 Posição da onda de choque no difusor ........................................................................... 92 5.4 Análise do Anteprojeto de um TVS ............................................................................... 95 5.5 Estimativa de onda de choque na seção de testes ........................................................... 100 6 CONCLUSÃO ......................................................................................................... 107 6.1 Conclusões e Sugestões .................................................................................................. 107 6.2 Proposta para Continuação da Pesquisa .................................................................... 110 REFERÊNCIAS .............................................................................................................. 111 APÊNDICE A .................................................................................................................. 115 18 1. INTRODUÇÃO 1.1 Definição do Problema Desde o fim da Segunda Guerra Mundial - quando foram criados os primeiros mísseis - até atualmente, o problema de controle desses continua sendo objeto de pesquisa. Essa observação não decorre somente do fato de as técnicas de defesa contra essas armas estarem cada vez mais apuradas mas também pela característica intrínseca de não linearidade e variância no tempo da dinâmica desses veículos. Essas não linearidades estão associadas, principalmente, aos regimes de velocidade (subsônico, transônico e supersônico) e atitude (arfagem, guinada e rolamento) de voo definidos para o envelope de operação desses veículos. Durante o desenvolvimento da tecnologia - foguete ou míssil - é prática comum estabelecer-se, inicialmente, parâmetros de referência de projeto que qualificam o veículo em um contexto denominado de “Ideal” (configuração “baseline”). Todos os refinamentos posteriores em termos de projeto de componentes (ou subsistemas) provocarão perdas nas características de performance/controle (por exemplo) quando comparadas com o modelo “Ideal”. Neste contexto, a elaboração de modelos acoplados preliminares, representativos de um dado veículo, é etapa fundamental para a linha de desenvolvimento de projeto que se pretende adotar desde que estes estabeleçam limites em termos do que se consegue alcançar. Adicionalmente, uma vez definido o projeto conceptual do veículo, torna-se possível estabelecer diretrizes na especificação de ensaios em túnel de vento, em voo ou ensaios ambientais, tendo em vista que as variáveis mais influentes estão definidas na modelagem matemática de cada subsistema. Dentro deste contexto, o Grupo de Pesquisa “Projeto Conceitual de Foguetes & Mísseis” do Instituto de Aeronáutica e Espaço (IAE) do Departamento de Ciências e Tecnologia Aeroespacial (DCTA) definiu um conjunto de atividades a longo prazo que 19 visam, entre outros aspectos, disponibilizar ferramentas de trabalho automatizadas para o auxílio no desenvolvimento de cada subsistema do veículo. Nesse conjunto de atividades a serem perseguidas, insere-se o desenvolvimento de modelos aerodinâmicos aplicados a configurações de veículos aeroespaciais. De forma geral, esse tipo de modelo representa o conjunto de coeficientes aerodinâmicos associados a cada ponto do envelope de operação (altitude, velocidade, atitude) do veículo. A determinação do modelo aerodinâmico de um veículo aeroespacial pode ser realizada de diferentes maneiras. Entre essas, destacam-se: métodos semi-empíricos (por exemplo, software Missile Datcom®), CFD - (“Computational Fluid Dynamics”), testes em voo e testes em túneis de vento. Ressalta-se, porém, que a metodologia adotada na definição final dos coeficientes aerodinâmicos é a de testes em túnel de vento. A primeira etapa para a realização desses testes em configurações de uso prático é a definição da matriz de testes (condições a serem ensaiadas) e, obviamente, a verificação da viabilidade de execução desta, tendo em vista que os custos associados neste empreendimento são elevadíssimos. Dentro desse enfoque, este trabalho disponibiliza um software de análise numérica de túneis de vento de Sopro (TVS), baseado na formulação matemática de “parâmetros concentrados”, onde as variações espaciais são desprezadas: propriedades (estados) do sistema são considerados homogêneos em todo o volume de controle, que tem como saídas principais a especificação geométrica dos respectivos módulos funcionais e o envelope de operações do TVS associados a um conjunto de requisitos operacionais fornecidos pelo usuário. Ressalta-se que, diferentemente do que ocorre com temas da área de aerodinâmica subsônica, existem disponíveis poucos dados confiáveis de geometrias e condições de funcionamento de túneis TVS tendo em vista o nível de sigilo envolvido na execução de testes em modelos (geralmente de aplicação militar); daí a grande aplicabilidade desta proposta de desenvolvimento. 20 1.2 Objetivo O presente trabalho tem por objetivos: � Anteprojeto aerotermodinâmico de um túnel de vento de sopro; � Desenvolvimento de uma ferramenta computacional para análise aerotermodinâmica e de controle de túneis de vento de sopro. 1.3 Túnel de Vento de Sopro Túnel de vento de Sopro (TVS) ou túnel de vento supersônico de sopro ou túnel de vento supersônico de disparo é uma instalação usada para testes aerodinâmicos em modelos de veículos aeroespaciais. O regime de escoamento do TVS desta tese é tipicamente supersônico. De um modo geral, o túnel TVS é caracterizado pelos seguintes módulos funcionais: (i) Reservatório de ar comprimido, (ii) Câmara de Tranquilização, (iii) Seção de testes e, finalmente, (iv) Difusor. Ressalta-se que se trata de uma descrição simplificada visto que também são utilizados em túneis TVS reais, compressor e secador de ar antes da entrada do reservatório de ar comprimido além de outros subsistemas associados ao controle do túnel, Figura 1.1. Os parágrafos subsequentes descrevem a funcionalidade de cada um dos módulos supracitados. (i) Reservatório de Ar Comprimido: o suprimento de ar principal de um TVS provém do reservatório de ar comprimido (Tanque de ar). Em algumas aplicações, esse módulo contém um regenerador de calor cuja finalidade principal é minimizar a queda da temperatura do ar durante o funcionamento do túnel. Do reservatório de ar comprimido, o ar escoa através de tubulação até uma válvula de bloqueio tipo gaveta. Junto à válvula de bloqueio, é instalada uma válvula de controle automático cuja finalidade é manter, durante 21 o escoamento, a pressão total constante na câmara de tranquilização ou o número de Reynolds constante na seção de testes através do controle de vazão de ar a partir do reservatório de ar comprimido. A seguir, o fluxo de ar passa por um difusor de grande ângulo cuja finalidade é reduzir a velocidade média do escoamento e recuperar parte da pressão dinâmica. (ii) Câmara de Tranquilização: o ar proveniente do módulo reservatório de ar comprimido passa por uma câmara de tranquilização (“Settling Chamber”) cuja finalidade principal é a redução da turbulência imposta ao escoamento durante a evacuação do reservatório de ar comprimido. Acoplada à câmara de tranquilização tem-se um bocal convergente-divergente (CD), algumas vezes com garganta ajustável e também chamado de “primeira garganta”, cuja finalidade é acelerar o escoamento até a velocidade supersônica desejada na seção de teste. Essencialmente, a relação de áreas da seção transversal da garganta e seção de testes, e as condições de estagnação na câmara de tranquilização definem o número de Mach na seção de testes. (iii) Seção de Testes: a seção de testes (ST) vem logo após o módulo câmara de tranquilização. É na seção de testes que se realizam os ensaios com modelos de aeronaves, mísseis, foguetes e configurações aerodinâmicas isoladas (como, por exemplo, asas, fuselagens e dispositivos aeronáuticos). (iv) Difusor: a seguir, o escoamento passa por um difusor CD que pode apresentar, também, garganta ajustável. A finalidade desse difusor é desacelerar o escoamento para as condições ambientes. A finalidade da garganta ajustável (“segunda garganta”) é auxiliar na estabilização da posição da onda de choque. Logo após o difusor está instalado um silenciador com descarga livre para a atmosfera. 22 Figura 1. 1- Desenho esquemático de um túnel de vento de sopro 1.4 Revisão bibliográfica A maioria das instalações de túneis de vento supersônico no mundo utilizam os mesmos procedimentos operacionais adotados nos anos 60 (MATSUMOTO, 2000). Os desenvolvimentos em andamento preocupam-se mais com particularidades e detalhes de ajuste fino para ganhos incrementais em eficiência e precisão do ensaio, tendo em vista o maior nível de exigência em termos de requisitos de projeto e os grandes custos associados à execução dos ensaios. Entre as referências mais importantes pesquisadas para a execução deste trabalho, destacam-se: Fung (1987) – Neste trabalho é apresentado o desenvolvimento e a implementação de um algoritmo de controle da pressão de estagnação do túnel de vento supersônico da Universidade Estadual da Pensilvânia. As maiores dificuldades na implementação deste sistema de controle foram as não linearidades inerentes à fenomenologia esperada para o campo de escoamento e as características da válvula de controle. Foi escolhido, para essa finalidade, um controlador PI de uma entrada e uma saída de sinal devido à sua simplicidade e disponibilidade. O desempenho do túnel de vento supersônico foi considerado satisfatório no sentido de que a pressão de estagnação foi controlada dentro de uma precisão aceitável em termos operacionais. Matsumoto (2000) – Foi desenvolvido um túnel de vento de Sopro para experiências em aeropropulsão no centro de pesquisas em aerodinâmica da 23 Universidade do Texas, em Arlington. O túnel é capaz de atingir regimes de velocidade variando de Mach 1.5 a Mach 4. O tempo útil máximo de corrida é da ordem de 2 segundos. A geometria da seção de teste é de 0,15 x 0,15 m (seção transversal). Novos componentes do túnel foram projetados e implementados com os componentes já existentes: bocal supersônico, tanque tubular de ar comprimido e compressor de alta pressão. Foi também desenvolvido um controlador pré-programado para sobrepujar a resposta lenta do controlador da válvula. Foi desenvolvido um perfil ideal de abertura da válvula para cada corrida com base nas condições específicas dos resultados dos ensaios e pressões iniciais e finais no tanque de ar comprimido. Após vários ensaios e respectivas interpolações corretivas, as perturbações de pressão na câmara de Tranquilização foram reduzidas para o valor típico de um por cento da pressão de estagnação. Foi feito um ensaio preliminar a Mach 2,5 e número de Reynolds de 63x106/m com uma pressão de estagnação de 928 kPa, e os resultados foram satisfatórios. Liu e Daley (2001) – São apresentados, neste trabalho, três esquemas de projeto de ajuste ótimo de controladores Proporcional-Integral-Derivativo (PID) para sistemas de controle industriais. Eles são controle PID de ajuste ótimo no domínio do tempo, controle PID de ajuste ótimo no domínio da frequência e controle PID de ajuste ótimo multi-objetivo. Esses esquemas podem prover parâmetros PID de ajuste ótimo de tal forma que as especificações desejadas para o sistema sejam satisfeitas mesmo no caso em que a dinâmica dele seja variável no tempo, ou os seus pontos de operação mudem. Eles são aplicados a três sistemas industriais: sistema hidráulico de controle de posição, sistema hidráulico de controle da rotação, e gaseificador. Bhoi e Suryanarayana (2008) – No túnel de vento de sopro NAL (National Aerospace Laboratories) 0,6 m, foi proposta a incorporação de um bocal de geometria variável através de um servo-atuador eletro-hidráulico para número de Mach variável de forma que o túnel possa partir de um baixo número de Mach e acelerar continuamente, até no máximo, Mach 4. Antes de finalizar a corrida no túnel, o servo- atuador eletro-hidráulico restabelece a configuração geométrica inicial do túnel. Trata-se de um trabalho puramente experimental. 24 Varghese e Binu (2009) – Foi desenvolvido um controlador Proporcional Integral (PI) para o controle efetivo da pressão de estagnação na câmara de Tranquilização de túnel de vento hipersônico de Sopro. Inicialmente, o projeto foi feito determinando-se os ganhos através de algoritmos de otimização (controle ótimo). A seguir, foram aplicadas técnicas de lógica difusa para melhorar a robustez e desempenho do controlador. Embora a formulação adotada tenha sido aplicada a um túnel de vento hipersônico, as equações descritivas são extremamente semelhantes às de um túnel de vento supersônico de Sopro. No cenário internacional da capacidade em ensaios em túnel de vento, o Brasil encontra-se em uma posição bastante desalentadora, uma vez que sua única instalação industrial é o túnel subsônico TA-2 do DCTA (FALCÃO,1996). Esse túnel de vento foi instalado na década de 50, chegando a desenvolver velocidades de até número de Mach 0,3. A seção de testes apresenta uma área útil de 2,0 x 3,0 m2. Esta realidade, ao longo dos anos de desenvolvimento de projetos na área aeroespacial, tem obrigado o Brasil a ensaiar seus modelos em túneis no exterior de modo limitado, seja por falta de recursos, seja por dificuldades oriundas do distanciamento entre a equipe de projeto e a equipe de ensaios em túneis. Há cerca de dez anos, o Ministério da Aeronáutica desenvolve esforços no sentido de projetar e construir um complexo de túneis transônico e supersônico de grande porte (Projeto TTS) dentro das dependências do Instituto de Aeronáutica e Espaço (IAE). O objetivo desse complexo é dar suporte ao atual estágio de desenvolvimento da indústria aeronáutica brasileira. Dentro deste contexto, destaca-se a construção do Túnel Transônico Piloto (TTP) em cujas instalações estão sendo envidados esforços no sentido de ser possível a realização de ensaios no baixo supersônico. Entretanto, instalações desta natureza envolvem grandes somas em termos de recursos humanos, financeiros e tecnológicos, o que não elimina a necessidade do Brasil utilizar as instalações de outros países. Entre os trabalhos desenvolvidos e publicados por pesquisadores brasileiros na área de TVS, destacam-se: Paglione (1978) – Nesta referência é apresentada uma interessante introdução à teoria do túnel aerodinâmico de sopro. Essa descrição engloba os princípios de funcionamento, formulação matemática que descreve o escoamento nos diversos 25 elementos do túnel e faz uma análise comparativa, baseada em informações bibliográficas, dos problemas associados à qualidade do escoamento durante o tempo de corrida, controle da pressão total do escoamento, tanto na partida como em regime quase permanente e determinação dos requisitos de precisão do sistema de controle. Essa dissertação de mestrado serviu como base para o projeto do túnel de sopro do ITA, em nível acadêmico. Silva, Falcão e Mello (2006) – Neste artigo, é desenvolvido um modelo matemático simplificado representativo de um TVS. Esse modelo considera o uso de regeneradores para controle de temperatura de estagnação do túnel de vento e válvula de abertura rápida para o controle de pressão de estagnação na câmara de tranquilização do túnel de vento. Constatou-se que o uso de um controlador PI para túneis, com dimensões semelhantes àquele apresentado em Fung (1987), conduziu a resultados de desempenho satisfatórios. Silva, Gamarra e Koldaev (2009) – O enfoque deste artigo foi o projeto de um sistema de controle PI para TVS e análise do subsistema difusor. São definidos dois controladores. O primeiro baseia-se na pressão de estagnação na câmara de tranquilização. O segundo é projetado com base na relação entre a pressão de estagnação e a temperatura na câmara de Tranquilização com finalidades de se controlar o número de Reynolds na seção de testes. O modelo matemático foi implementado em blocos Simulink®. O desempenho do túnel de vento supersônico foi considerado satisfatório, em comparação com testes reais, com túnel de vento de sopro, feitos na Pennsylvania State University em 1987. Uma equipe de especialistas do DCTA, da Força Aérea Brasileira (FAB), está agora na África do Sul participando do desenvolvimento conjunto do míssil ar-ar A- Darter (AAM). O contrato é entre o Ministério da Defesa do Brasil e a FAB, por um lado e, de outro, a Agência de Desenvolvimento e Compra de Materiais de Defesa, ARMSCOR, organização ligada ao Departamento de Defesa da África do Sul. Trata-se de um acordo de alto nível entre os respectivos governos. Os especialistas brasileiros 26 provenientes da FAB, empresa MECTRON, empresa AVIBRÁS e OPTOELETRÔNICA estão trabalhando com seus parceiros sul-africanos, locados nas dependências da empresa DENEL Aerospace, no desenvolvimento desse míssil. O processo de capacitação proposto visa estabelecer um mecanismo de comunicação que permita a Transferência de Tecnologia que pode ser aproveitada de forma direta pelo setor Aeroespacial do IAE. Entre os conhecimentos adquiridos na linha de pesquisa aerodinâmica, destacam-se a capacitação no dimensionamento de matriz de testes e a realização de ensaios em diferentes túneis transônicos e supersônicos (ONERA - França, CSIR - África do Sul). Dessa experiência, pôde-se constatar que o entendimento preliminar de todo o procedimento de teste e estimativa de resultados para um dado modelo de ensaio resulta em excepcionais reduções de custo para o projeto. Esta tese se enquadra diretamente nessa linha de pesquisa e desenvolvimento. 1.5 Escopo do Trabalho O Capítulo 2 apresenta a formulação matemática aplicável a cada um dos subsistemas do TVS. São definidas as formulações teóricas conjuntamente com as respectivas hipóteses adotadas e os principais aspectos teóricos do controlador PI, projetado para o TVS. Neste mesmo capítulo é disponibilizado a metodologia numérica de integração das equações diferenciais ordinárias (EDO), não lineares, variantes no tempo, que caracterizam a planta, neste caso, o TVS. São apresentadas também as principais características do código computacional desenvolvido na presente abordagem. O Capítulo 3 detalha o módulo “Seção de Testes” de um TVS. São apresentadas a metodologia desenvolvida para a estimativa da forma e localização da onda de choque que se forma em um modelo colocado no interior da seção de testes do TVS virtual e a interface desenvolvida com o código MISSILE DATCOM para o cumprimento desta estimativa. No Capítulo 4 é apresentada a metodologia de dimensionamento de um TVS que será utilizado nas simulações deste trabalho. Também são apresentados os principais 27 limites no uso da formulação adotada quando da aplicação em problemas de interesse prático. No Capítulo 5 é mostrada uma seleção representativa de resultados que dizem respeito à validação da metodologia, e resultados de interesse geral. Finalmente são apresentados, no capítulo 6, comentários referentes à aplicabilidade do presente desenvolvimento, às suas limitações e são discutidas algumas recomendações para trabalhos futuros. 28 2. FORMULAÇÃO MATEMÁTICA PARA O TVS Neste tópico é apresentada a formulação matemática associada ao funcionamento de um TVS em regime quase estacionário e a metodologia numérica de resolução das equações diferenciais. A Figura 2.1, a seguir, apresenta as variáveis de estado a serem empregadas na formulação dinâmica. Figura 2. 1 - Desenho esquemático de um túnel de vento de sopro 2.1 Hipóteses A análise dinâmica de TVS é dividida nos módulos, quais sejam: reservatório de ar comprimido, câmara de tranquilização, bocal, seção de testes e difusor. Na metodologia definida para este trabalho, cada subsistema é representado por volumes de controle (parâmetros concentrados) cuja formulação matemática é baseada nas equações de conservação de massa e energia. Ressalta-se que esta aproximação representa cada subsistema do túnel de vento como uma “caixa preta”, ou seja, esta metodologia não contempla a presença de turbulência e efeitos locais de fricção e transferência de calor no campo de escoamento no interior destes volumes de controle, 29 e sim, apenas as influências globais. Assume-se que a distribuição de temperaturas, pressões e massa específicas seja uniforme em cada face do volume de controle representativo dos subsistemas da planta em estudo. O uso da técnica de parâmetros concentrados é justificada em função da praticidade de obtenção de resultados. Essa praticidade é extremamente útil no estudo e seleção de sistemas de controle a serem desenvolvidos para a planta em estudo bem como na determinação de parâmetros globais de desempenho, tais como, tempo de execução de testes, variações de pressão e temperatura ao longo de cada volume de controle, entre outros aspectos. A variação de energia potencial do gás também é desprezada. Considera-se que o fluido de trabalho comporta-se como um gás, térmica e caloricamente perfeito. Finalmente, é importante ressaltar que a metodologia proposta não leva em consideração a inicialização (“start up”) do túnel, desde que o principal objetivo é o estudo das características de desempenho do túnel de vento. 2.2 Reservatório de Ar Comprimido O reservatório de ar comprimido é o primeiro módulo da planta a ser considerado para o balanço de massa e energia. É adotada a hipótese de que a contribuição de massa do compressor seja desprezível durante a corrida do túnel. Portanto, a taxa de diminuição de massa do tanque é igual à taxa de saída de massa através da válvula de controle, ou seja: v T T m Vdt d � 1 �� � , (2.1) sendo T� a massa específica do ar no reservatório de ar comprimido, vm� a massa de ar que passa através da válvula de controle e TV o volume do tanque de ar comprimido. A variável t representa o tempo. O subscrito “T” refere-se ao reservatório de ar comprimido. 30 Assumindo que a perda de energia através da válvula seja desprezível, pode-se arguir que a variação de energia interna( TU ) no reservatório de ar comprimido é dada por: 2 2 1 vvvv T vmhm dt dU �� ��� . (2.2) Nesta formulação vh representa a entalpia específica do ar que passa através da válvula de controle e vv é a velocidade do ar nesta válvula. Pode-se escrever a equação (2.2) em termos das condições de pressão( TP ) e temperatura( TT ) de estagnação no reservatório de ar comprimido (FUNG, 1987): v T TT m V RT dt dP ��� � � �� � � �� � . (2.3) O parâmetro � é a razão de calores específicos ( vp cc ) e R é a constante de gás. É instrutivo ressaltar que a válvula de controle utilizada na planta é o modelo Camflex DN 8”, da Fisher Controls Company (1984). A principal razão de usá-la como referência, neste trabalho, é a disponibilidade de informações no trabalho de Fung (1987). O fluxo de massa definido para diferentes posições do obturador da válvula de controle é calculado através da relação empírica (2.4). �� � � �� � � � � � T Tg T v P P senPC T m 71.2295810.2 8 � . (2.4) O parâmetro gC é função da posição do obturador da válvula de controle, ou seja, � ��gg CC � , sendo � o ângulo de abertura da válvula. P� é a diferença entre as pressões de estagnação no reservatório de ar comprimido( TP ) e na câmara de tranquilização ( 0P ). A Tabela 2.1 indica a relação funcional � ��gg CC � . 31 2.3 Câmara de Tranquilização O segundo volume de controle de interesse é aquele representativo da câmara de tranquilização. O ar comprimido, proveniente do reservatório, entra na câmara de tranquilização através da válvula de controle, passa para o bocal CD e, daí, para a seção de testes. Tabela 2. 1– Relação entre gC e � � [graus] gC 0 0 10 194 20 1680 30 3767 40 6230 50 9288 60 12835 70 16351 80 18942 90 23120 Desta forma, um balanço da massa para esse segundo volume de controle fornece a equação (2.5): � �tv mm Vdt d �� �� 0 0 1� . (2.5) O subscrito “0” se refere à câmara de tranquilização e, o subscrito “t” (“throat”), às condições de escoamento saturado na garganta do bocal de entrada da seção de testes. A energia que entra no volume da câmara de tranquilização com a vazão de massa vm� menos a energia que sai através do bocal com a vazão de massa tm� é igual à taxa de variação de energia( 0U ) dentro da câmara de tranquilização, ou seja: 32 220 2 1 2 1 ttttvvvv vmhmvmhm dt dU ���� ���� . (2.6) Analogamente ao procedimento adotado para o reservatório de ar comprimido, pode-se escrever a equação (2.6) em termos das condições de estagnação, qual seja: � �0 0 0 TmTm V R c c dt dP tTv v p �� ��� � � �� � � � . (2.7) 2.4 Bocal do TVS O bocal do TVS em estudo é simétrico, convergente-divergente (CD) e, conforme dito anteriormente, o escoamento é considerado isentrópico. Desde que essa região é a menor área de passagem de toda a planta, considera-se que ela permaneça saturada durante toda a simulação. Considerando-se que o ar se comporta como um gás perfeito e que o estado de referência é o estado de estagnação, pode-se escrever a equação da vazão mássica durante a ocorrência de escoamento sônico na garganta do bocal como sendo: Dtt CAPm 0�� , (2.8) sendo DC o coeficiente de descarga do bocal e tA a área da seção transversal na região da garganta do bocal. A equação (2.9) exibe a expressão de DC . � �12 1 2 1 1 2 � � �� � � �� � � ��� � � �� � � � � � � � T D RT C . (2.9) ‘ A partir dessas relações, é possível demonstrar que a definição do número de Mach na seção de testes ( STM ) é, essencialmente, dependente apenas da relação de áreas da seção transversal da seção de testes ( 0A ) e da primeira garganta do túnel ( tA ), desde que o escoamento esteja saturado na garganta do bocal de entrada, ou seja: 33 � � � �� �12 1 2 2 1 2 11 � �� � � � � � � ! " � � � � � � � � � � � � � � � � � � � � ST ST ST t M M A A . (2.10) 2.5 Seção de Testes e Difusor Não foi implementada uma formulação específica para a seção de testes uma vez que os parâmetros que a caracterizam são obtidos indiretamente a partir da solução dos outros subsistemas. Entretanto, a definição da pressão de estagnação na seção de testes ( STP ) como condição de projeto deve ser enfatizada. É através desta condição que será posicionada a onda de choque no difusor e, em consequência, as condições operacionais do túnel de vento. Uma vez que um TVS esteja em operação, é necessário que o escoamento proveniente da seção de testes seja devolvido nas condições ambientes na saída do túnel de vento. O projeto da maioria dos túneis de vento que trabalham com escoamentos supersônicos inclui, para esta finalidade, um dispositivo denominado de difusor. Este subsistema é composto de um bocal, cuja menor seção transversal é denominada de “segunda garganta”. A eficiência de um difusor está relacionada diretamente à sua geometria desde que para uma dada condição funcionamento ter-se- á uma pressão recuperada como função do trem de ondas de choque oblíqua e da intensidade da onda de choque normal que ocorre em seu interior, Figura 2.2. O número de Mach em que ocorre a onda de choque normal no difusor deve ser menor que aquele dimensionado para a seção de testes. A condição de projeto para a alimentação de um túnel é baseada na maior perda de carga que ocorre em seu interior, a qual, para o caso de túneis supersônicos, deve ocorrer na seção de testes, POPE e GOIN (1965). É instrutivo ressaltar também que a estimativa do tempo de ensaio é fundamentalmente dependente da posição da onda de choque no difusor desde que, sob certas condições de relação de pressão, o escoamento no interior do difusor fica instável e, em consequência, esta onda de choque pode ser “engolida” pelo difusor e comprometer as condições estabelecidas na seção de testes. 34 (a) Região em que ocorre a onda de choque em tubo de área constante (M=2,06) (b) Detalhamento da região em que ocorrem as ondas de choque (M=2,55) Figura 2. 2 - Fotografia Schlieren de um difusor supersônico (Neumann, Lustwerk,1949) Neste trabalho é considerado que o bocal do difusor seja convergente-divergente. Dentro deste enfoque, o projeto do difusor é baseado na área da segunda garganta que produz uma onda de choque de intensidade SdM na porção divergente, Figura 2.2 conforme (NEUMANN, LUSTWERK, 1949). Neste tópico é apresentada a formulação matemática implementada para um difusor ideal. Não obstante, ressalta-se que as condições estabelecidas na seção de testes são calculadas com base em dados experimentais de forma a garantir o estabelecimento de uma onda de choque o mais próximo possível da garganta do difusor. Dentro deste enfoque, o problema pode ser colocado conforme segue: Dados que devem ser disponibilizados para o projeto do difusor: i. � : Perda em pressão de estagnação entre o bocal de aceleração (1ª garganta) e a entrada do difusor (2ª garganta); 35 ii. td ed A A : Razão de Áreas: Área de saída do difusor por Área da 2ª garganta; iii. dp : Pressão estática na saída do difusor. Parâmetros calculados a partir dos requisitos supracitados: i. tdA : Área da garganta do difusor (2ª garganta); ii. td s A A : Posição da onda de choque logo após a segunda garganta (caso Ideal); iii. SdP : Pressão de estagnação antes da onda de choque no difusor. A Figura 2.3 abaixo apresenta exemplos de difusor com geometria fixa. Figura 2. 3 – Exemplos de Difusor Supersônico (Ref.: Neumann, Lustwerk, 1949) 36 A área da 2ª garganta pode ser calculada a partir do princípio da continuidade, qual seja: DdtdSTDtt CAPCAPm �� 0� (2.11) Considerando que a temperatura de estagnação não se altera significativamente durante o funcionamento do túnel e que existe uma perda de pressão de estagnação entre a 1ª garganta e a segunda garganta (� ), tem-se: ,1 0 00 t t ST t td A P AP P AP A �� ��� (2.12) note que 1#� . O número de Mach na saída de um difusor ideal pode se calculado conforme segue (JOHN e KEITH, 2006): 22 1 12 2 1 2 1 2 1 1 1 1 �� � � �� � � �� � � �� � � �� � � �� � � ��� � � �� � � � ��� � � �� � � � � � �� � � ed td ed Sd ed A A p P M � � ���� . (2.13) O subscrito “ed” se refere à saída do difusor e o subscrito “td” se refere à área da segunda garganta. Ressalta-se que todos os subscritos terminados em “d” se referem ao difusor. A variável SdP representa a pressão de estagnação antes da onda de choque no difusor. O procedimento adotado para se calcular a posição da onda de choque segue os seguintes passos: i. Determinar edM a partir da Eq.(2.13) considerando-se 0PPSd �� ; ii. Com edM determinar d ed P p , desde que 12 2 11 � � � � � � � � � �� � � � ed d ed M P p . Note que dP é a pressão de estagnação no difusor após a onda de choque; 37 iii. Desde que edM < 1, pode-se adotar atmed pp � . Em consequência, ed d Sd ed Sd d p P P p P P � ; iv. Conhecendo-se Sd d P P , determina-se o número de Mach antes da onda de choque( SdM ) através das relações de salto; v. Finalmente, td s A A pode ser calculado através de uma equação- com a mesma estrutura matemática da Eq.(2.10). Note que sA é a área da seção transversal do difusor onde a onda de choque está localizada e tdA é a área da seção transversal da segunda garganta. O procedimento supracitado foi implementado numericamente no desenvolvimento do programa computacional que simula a dinâmica de um TVS. 2.6 Sistema de Controle A razão de se instalar um controlador em TVS é a necessidade de se garantir a qualidade do escoamento (uniformidade dos números de Mach e Reynolds ) na seção de testes do túnel. Esta qualidade é avaliada pelo grau de uniformidade do escoamento e pelas propriedades termodinâmicas definidas, a priori, nos requisitos de ensaio. Um exemplo típico de resultados alcançados em túneis controlados por válvulas reguladoras está citado na referência, devido a MARVIN (1987). É possível se estimar os coeficientes de arrasto e de pressão em modelos de aeronaves, com erros inferiores a 1% relativos à condição de projeto. Para atingir estes critérios, faz-se necessário garantir uniformidade do número de Mach na seção de testes de $ 2,2 % (neste exemplo, Mach = 3). Da mesma forma é difícil manter o número de Reynolds constante durante a corrida, visto que a temperatura do tanque de ar comprimido cai continuamente. Requisitos semelhantes ao exemplificado é parte fundamental do contrato que se estabelece entre a empresa desenvolvedora de um veículo aeronáutico (cliente) e a empresa que executará os ensaios. 38 2.7 Controlador O controlador projetado para este desenvolvimento é do tipo Proporcional- Integral (PI). Esta opção foi baseada nas diversas referências consultadas na área e mencionadas no Capítulo 1 desta tese. Este controlador possui dois parâmetros livres, quais sejam: ganho proporcional ( pK ) e ganho integral ( iK ). Estes ganhos devem ser estimados de forma a garantir a estabilidade de operação e minimizar a duração do transiente inicial a fim de se obter uma corrida estável e a mais longa possível. Essencialmente, o sistema de controle trabalha da seguinte maneira: considere que o controlador se baseia na pressão de estagnação da câmara de tranquilização; o controlador digital da válvula de controle compara o valor real da pressão de estagnação com o valor definido em requisito (valor de referência) e envia um sinal corretivo de vazão mássica de ar que é função dos ganhos pK e iK . A função de transferência de um controlador PI é dada pela equação (2.14), OGATA (1997). � � � � � � �� � � �� � � ��� sK K sE s sG i p 11� , (2.14) sendo � �s� a posição de abertura da válvula de controle e � �sE o erro de sinal entre a entrada de referência e a saída da planta. Neste trabalho estão sendo apresentados 2 (dois) sistemas de controle para TVS. O primeiro tem como ponto de projeto (“set point”) a pressão de estagnação na câmara de tranquilização e o segundo utiliza o número de Reynolds, definido em requisito para a seção de testes, como ponto de projeto. Ambos controladores estão descritos em maiores detalhes nas seções subsequentes. 39 2.7.1 Controlador baseado na pressão de estagnação da câmara de Tranquilização Para este tipo de controlador, � �s� é a posição de abertura da válvula de controle e � � � � � �sP sP sE ojetoPr 0 01�� representa o erro de sinal entre a entrada de referência ( 1int 0 �setpoP ) e a saída da planta � � � �sP sP ojetoPr 0 0 que representa a pressão real medida. Note que o parâmetro � �sP ojetoPr 0 representa a pressão de estagnação definida como ponto de projeto na câmara de tranquilização e � �sP0 é a pressão de estagnação variante no tempo na câmara de tranquilização. Aplicando-se a transformada inversa de Laplace entre a entrada e a saída � �t� do controlador PI, pode-se escrever a equação (2.14) como sendo: � � � � � � � � � �� � � � �� � � �� �� � � �� � � �� tP tP K K dt tP tP d K dt td ojeto i p ojeto p Pr 0 0 Pr 0 0 1� , (2.15) 2.7.2 Controlador baseado no número de Reynolds da seção de testes Com base na discussão precedente, pode-se dizer que é possível controlar as condições na seção de testes, utilizando-se como ponto de projeto a pressão de estagnação na câmara de tranquilização. Este é o procedimento normalmente adotado por inúmeras instalações industriais, FUNG (1987). Contudo, é de se esperar uma significante queda na temperatura de estagnação em toda a planta durante o processo de esvaziamento do reservatório de ar. Dependendo da vazão mássica na planta e do volume do reservatório, esta variação em temperatura altera significativamente o número de Reynolds ( Re ) configurado para o ensaio. Neste contexto, um sistema de controle baseado no número de Reynolds na seção de testes deve ser analisado de forma a se avaliar a possibilidade de utilizá-lo como parâmetro de referência. Este trabalho propõe uma metodologia de controle que permite minimizar o efeito da queda de temperatura no número de Reynolds estabelecido para a seção de testes. Para tanto, é necessária uma tomada de pressão 40 (localizada na câmara de tranquilização) e uma tomada de temperatura estática (localizada na seção de testes). Os parágrafos subsequentes descrevem esta metodologia. A partir da definição de um processo isentrópico, fazendo-se STPP �0 , pode-se afirmar que: 112 2 11 �� �� � � � � � � �� � � � � � FM p P ST ST ST (2.16) Ressalta-se que na Eq. (2.16), STP é a pressão total e STp é a pressão estática na seção de testes. Sendo ST STT F % � e ST% a temperatura estática na seção de testes. Desde que se está adotando a hipótese de gás termica e caloricamente perfeito, tem-se: 1 1 � � � � FRT P ST ST ST . (2.17) A partir da definição do número de Reynolds, pode-se deduzir que: � � ST STSTSTST RTFMD t � �� � � 12 1 )Re( � � � (2.18) sendo ST� a viscosidade dinâmica do gás na seção de testes e STD o comprimento de referência adotado no cômputo do número de Reynolds. Neste caso, STD foi adotado como sendo o diâmetro da seção de testes. A partir das definições: 4 2 ST ST D A & � (2.19) e 41 ST STST STST RT P F P RT � �� � 1 1 � � , (2.20) pode-se escrever o número de Reynolds como função das condições de estagnação na seção de testes, qual seja: ST ST T P t 5.1 )Re( � , (2.21) sendo F 1 2 � , 08 1 103110.2 �� x e R MA STST & � 2 2 4 � . Analogamente ao desenvolvimento apresentado no item 2.7.1, obtém-se a equação do controlador: � � � � � � � � � � � � �� � � � � � � �� ojeto i p ojeto p t K K dt t d K dt td Pr Pr Re Re1Re Re � . (2.22) 2.7.3 Determinação dos ganhos do controlador Para o caso em estudo, não existe uma fórmula específica e direta para se calcular os ganhos associados ao controlador. Essencialmente, os ganhos são determinados através de metodologias já consagradas na literatura. A metodologia utilizada aqui é baseada na referência devido a Varghese e Binu (2009). A estrutura matemática das equações usadas para a determinação dos ganhos da planta daquele artigo é muito semelhante às equações do TVS deste trabalho. A determinação dos ganhos Kp e Ki , conforme apresentado no Apêndice A, envolve o cumprimento das seguintes etapas: a. Linearização das equações (2.3), e (2.7). Utilizou-se a metodologia adotada por CHAU (2001) neste tópico; 42 b. Construção das matrizes das equações de espaço de estados; c. Determinação das funções de transferência; d. Determinação dos ganhos Kp e Ki do controlador; Nos itens subsequentes são apresentados os desenvolvimentos supracitados: � Linearização das equações (2.3) e (2.7): Utilizou-se a metodologia adotada por CHAU, P. C., (2001) neste tópico; a Tabela 2.2 sintetiza o que foi feito. Os parâmetros k1 à k5 são constantes decorrentes do processo de linearização. O superscrito “L” nas expressões linearizadas é utilizado apenas para informar que se trata de uma expressão linearizada. A variável L g� é o ângulo de abertura da válvula de controle. Tabela 2. 2– Equações linearizadas: Reservatório de ar e Câmara de tranquilização Eq. Expressão não linearizada Expressão linearizada Eq. (2.3) v T TT m V RT dt dP ��� � � �� � � �� � L t L g L t Pkk dt dP *21 �'� � (2.23) (2.7) � �0 0 0 TmTm V R c c dt dP TTv v p �� ��� � � �� � � � LL t L g L PkPkk dt dP 0543 0 * '��'� � (2.24) � Construção das matrizes das equações de espaço de estados: Para este tipo de problema é pratica comum se utilizar a formulação relativa a modelos no espaço de estados. Para o uso desta formulação é necessário colocar o sistema de equações diferenciais da planta (túnel de vento) na seguinte forma: BuAxx ��� , e (2.25) DuCxy �� (2.26) O diagrama de blocos correspondente às equações (2.25) e (2.26) está ilustrado na figura seguinte. 43 Figura 2. 4 – Espaço de Estados - Diagrama de Blocos Na Figura 2.4 os valores iniciais estão identificados por x0. A variável x é o vetor das variáveis de estado, u é o vetor da excitação de entrada (no caso a variação do ângulo �), y é o vetor de saída. As matrizes A, B, C e D representam a matriz da planta matriz de entrada, matriz de saída e matriz de transmissão direta, respectivamente. De fato, poucos processos e sistemas possuem características tais que o sinal de entrada consiga agir diretamente sobre o sinal de saída. Devido à isto a matriz D é usualmente zero. No caso em estudo, a matriz D é também igual a zero porque não existe nenhum componente do sinal de entrada que tenha influência direta sobre o resultado de saída. A construção das matrizes no espaço de estados é feita da seguinte forma, adotando para as equações (2.23) e (2.24) : 1x dt dP L t �� ; 1xPL t � ; 2 0 x dt dP L �� ; 20 xP L � ; 1uL g �� (2.27) Aplicando as expressões (2.27) nas equações (2.23) e (2.24), tem-se: 12111 * xkukx �'�� (2.28) 2514132 * xkxkukx '��'�� (2.29) Pode-se reescrever (2.28) e (2.29) como: 112121 0.0* ukxxkx '�'��� (2.30) 44 1325142 * ukxkxkx '�'��� (2.31) ou seja: 1 3 1 2 1 54 2 2 1 0 u k k x x kk k x x � � � ! " � ( ) * + , - � � � ! " � ( ) * + , - � � (2.32) � � ( ) * + , - � ( ) * + , - 2 1 2 1 10 x x y y , (2.33) Lembrando que não está sendo considerado ruído, ou seja: 0�D (2.34) Comparando-se as equações (2.32), (2.33) e (2.34) com as equações (2.25) e (2.26), pode-se escrever: � � � ! " � 54 2 0 kk k A ; � � � ! " � 3 1 k k B ; � �10�C ; D = 0 (2.35) � Determinação das Funções de Transferência (FT): A função de transferência que corresponde ao sistema de equações (2.25) e (2.26) é dada por (VON ZUBEN, 2003): � � � � � � DBAIC ���� �1)( s sU sY sGSS (2.36) Na equação (2.36), I é a matriz identidade e as demais matrizes foram definidas anteriormente. Segundo XUE, CHEN e ATHERTON (2007), é possível determinar os polos e zeros da planta através do comando MATLAB®: A = [k2, 0; k4, k5]; B = [k1; k3]; C = [0, 1]; D = 0; G = ss(A,B,C,D);G1 = tf(G) A expressão da FT resultante para a formulação com polos e zeros é dada por: 45 � � � � � �21 1)(1 psps zs ksG ��� � �� (2.37) Nesta expressão, k é o ganho da planta, 1z é o zero, 1p é o primeiro polo e 2p é o segundo polo. Esta função de transferência pode ser melhorada a partir da inclusão de um atraso na resposta do sistema (VARGHESE e BINU, 2009), definido por � �sGD , e de um atuador para a válvula de controle (CHAU, 2001), representado por � �sGv . Com estas aproximações, a função de transferência global da planta é dada por: � � � �sGsGsGsPTF vD ''� )()(1 (2.38) sendo: � � 1 1 �� ��� � s s sGD % % , (2.39) e, � � 1� � s k sG v v v % (2.40) Na Eq.(2.40), vk é o ganho do atuador da válvula e v% é a constante de tempo do atuador da válvula. Para um controlador da pressão de estagnação, Eq.(2.14), tem-se: � � � � � � �� � � �� � � ��� sK K sE s sG i p 11� . (2.41) A expressão (2.41) pode também ser reapresentada como: � � � � � � s sK K K sK KsGsG i i p i pC 111 � '�� � � �� � � ��� � � �� � � ��� , (2.42) ou seja: 46 � � � � s sK K K sG i i p C 1� '�� � � �� � � � , (2.43) Sendo que as variáveis pK e iK são, respectivamente, os ganhos proporcional e integral do controlador, a serem determinados. � Determinação dos ganhos pK e iK do controlador: O diagrama de blocos em malha fechada com o controlador e a planta (Túnel de Vento) está apresentado na Figura 2.5: Figura 2. 5 - Diagrama de blocos em malha fechada para o TVS Para a determinação dos ganhos pK e iK do controlador para o TVS usou-se como valores de partida os ganhos do controlador usado no trabalho de VARGHESE e BINU (2009). Os referidos ganhos foram: VB pK = 2,2676e-007 e VB iK = 3,118562018. Ressalta-se que o superescrito “VB” indica que estes valores estão associados àqueles definidos por VARGHESE e BINU (2009). De acordo com a Figura 2.5: � � � � � � � � � �sPTFsG sPTFsG sG C C T '� ' � 1 (2.44) Na equação (2.44), � �sGT é a função de transferência, para o túnel completo incluindo o controlador e a realimentação (“feedback”). A metodologia adotada 47 consiste em excitar a função de transferência da planta (Eq. 2.44) com entrada degrau unitário e verificar se a resposta no tempo é estável e se as demais características tais como sobressinal máximo (“overshoot”), tempo de subida, tempo de pico e tempo de estabilização são coerentes para o caso específico. Esta coerência é avaliada com base em dados históricos. Os ganhos pK e iK do controlador são selecionados a partir da comparação entre a resposta da planta em estudo e valores de referência (conjuntamente com a experiência do grupo de trabalho). O MATLAB® dispõe de vários recursos que tornam esta atividade bastante prática. Entre estes recursos, destaca-se a ferramenta SISOTOOL. Mais detalhes de todas as etapas para a determinação dos ganhos do controlador, incluindo os respectivos valores numéricos, constam no Apêndice A. Seguindo a metodologia supracitada, os valores finais aceitáveis encontrados para os ganhos pK e iK do controlador foram: pK = 226,76e- 007 e iK = 0,3118562. Conforme exposto no Apêndice A, decidiu-se investigar também o comportamento do TVS com os ganhos do controlador usado por FUNG, 1987, onde: pK = 2,0 e iK = 0,08, e descobriu-se que, ao fazer-se a comparação de desempenho entre estes dois controladores, o controlador usado por FUNG, 1987, apresentou um desempenho muito melhor do que o controlador oriundo do trabalho de VARGHESE e BINU, 2009. Devido a isso decidiu-se usar nesta tese os ganhos do controlador FUNG, 1987. A comparação dos resultados usando os dois controladores é apresentada no capítulo 5. 2.8 Implementação Numérica A partir dos itens anteriores foi possível derivar um modelo matemático que represente algumas das principais características necessárias para o projeto conceitual de um túnel de vento de sopro. Esse modelo matemático está sintetizado na Tabela 2.3. As equações que representam o TVS são equações diferenciais ordinárias (EDO), não lineares e dependentes do tempo. Este sistema tem como variáveis de estado a pressão de estagnação e a massa específica do ar no reservatório de ar comprimido 48 ( TP , T� ), a pressão de estagnação e massa específica do ar na câmara de tranquilização ( 0P , 0� ) e o ângulo de abertura da válvula de controle (� ). Trata-se de um conjunto de 5(cinco) equações diferenciais e 2(duas) relações constitutivas utilizadas para fechar o sistema, quais sejam: vazão mássica que passa pela válvula de controle ( vm� ) e vazão mássica que passa através do bocal associado à 1ª garganta ( tm� ). Os parâmetros livres necessários para a resolução deste sistema, estão listados a seguir: i. Gás: razão de calores específicos e constante de gás (� , R ); ii. Válvula de Controle: relação funcional entre a constante da válvula ( gC ) e o ângulo de abertura da válvula (� ); iii. Sistema de Controle: ganhos proporcional ( pK ) e integral ( iK ). 2.9 Metodologia de integração Neste trabalho é utilizado a função (“solver”) ODE15s do aplicativo MATLAB® para a resolução do sistema de EDO representativo do TVS. Esta função é normalmente utilizado para a resolução de problemas com variações acentuadas das variáveis (“stiff problems“) e, para o caso em questão, é utilizada a opção de passo no tempo variável. O sistema de EDO foi implementado na forma de blocos utilizando-se o software SIMULINK®. A Tabela 2.4 exibe o descritivo da entrada de dados para o programa principal. A Figura 2.6 exibe o diagrama utilizado para o controlador de pressão de estagnação da câmara de tranquilização. Uma estrutura análoga é adotada para o controlador de número de Reynolds na seção de testes. Alem disso também foi feita uma segunda versão de programa em Matlab® sem o uso do SIMULINK®.. Ressalta-se que a verificação da implementação numérica foi realizada a partir do desenvolvimento de um segundo código computacional, escrito em MATLAB®, cuja finalidade principal foi detectar possíveis discrepâncias nos resultados obtidos com o código computacional escrito em SIMULINK®. A versão escrita em MATLAB®, 49 embora não apresentasse uma interface gráfica para entrada de dados, demonstrou ser mais rápida em termos de tempo de processamento. 50 Tabela 2. 3– Modelo matemático para o TVS MÓDULO MODELO MATEMÁTICO REFERÊNCIA Reservatório de Ar Comprimido v T TT m V RT dt dP ��� � � �� � � �� � v T T m Vdt d � 1 �� � (2.3) (2.1) Válvula de Controle � � � � � � � � � � � T Tg T v P P PC T m 71.2sin295810.2 8 � (2.4) Câmara de Tranquilização � �0 0 0 TmTm V R c c dt dP tTv v p �� ��� � � �� � � � � �tv mm Vdt d �� �� 0 0 1� (2.7) (2.5) Bocal (1ª garganta) Dtt CAPm 0�� (2.8) Sistema de Controle: Pressão de Estagnação � � � � � � �� � � �� � � �� �� � � �� � � �� ojeto i p ojeto p P tP K K dt P tP d K dt td Pr 0 0 Pr 0 0 1� (2.15) Sistema de Controle: Número de Reynolds � � � � � � � � � � � � �� � � � � � � �� ojeto i p ojeto p t K K dt t d K dt td Pr Pr Re Re1Re Re � (2.22) 51 Tabela 2. 4– Entrada de Dados MÓDULO PARÂMETRO OBSERVAÇÃO Reservatório de Ar Comprimido e Câmara de Tranquilização Volume Pressão e Temperatura Inicial Sistema de Unidade: SI Seção de Testes Geometria Pressão e Temperatura Inicial Condições de Projeto As condições de projeto são definidas pelo números de Mach e Reynolds Controlador Ganhos Proporcional e Integral Curva de vazão da válvula de controle Adimensional Difusor Perda em Pressão de Estagnação Razão de Áreas Pressão estática: saída do difusor � td ed A A dp Figura 2. 6 - Diagrama SIMULINK: Controlador de Pressão de Estagnação 52 3. SEÇÃO DE TESTES Com a implementação numérica do sistema de equações diferenciais ordinárias representativo do comportamento aerotermodinâmico da planta em estudo (Tabela 2.3) foi possível gerar o código computacional WINDT (atualmente na versão 03). Essencialmente, este código computacional inclui todos os módulos representativos de um túnel de vento de sopro convencional com exceção, da inclusão de um modelo na seção de testes. Obviamente, a implementação deste recurso exige uma formulação matemática mais robusta que aquela utilizada por este trabalho como, por exemplo, técnicas de Dinâmica dos Fluidos Computacional. A fim de viabilizar uma primeira análise, em termos de engenharia, procurou-se trabalhar no desenvolvimento de uma interface entre o programa WINDT e o código comercial MISSILE DATCOM®. Esta interface foi idealizada visando o cumprimento de três objetivos: (i) determinar os coeficientes aerodinâmicos associados a um modelo colocado no interior da seção de testes do túnel de vento projetado como se tivessem sido obtidos através de medições; (ii) estimar o posicionamento da onda de choque curva quando o modelo estiver a zero grau de ângulo de ataque e (iii) realizar um cálculo preliminar de cargas sobre o modelo. Neste trabalho foram implementados os itens (i) e (ii). A interface entre os referidos códigos é projetada para permitir o uso de dados da simulação do TVS na execução do código MISSILE DATCOM®. É importante ressaltar que se tratam de dois módulos distintos, quais sejam, o programa que calcula as condições aerotermodinâmicas de um túnel de sopro e um programa que calcula coeficientes aerodinâmicos para configurações do tipo míssil com se estivessem sendo lidos pela instrumentação do TVS. Os dados de saída de um programa (WINDT) servem como dados de entrada para o outro programa (MISSILE DATCOM®) e o processo de execução ocorre de forma automática. Não existe um acoplamento matemático entre as referidas formulações. A implementação supracitada permite que o usuário do programa WINDT obtenha de forma direta todos os dados aerodinâmicos associados à simulação de um modelo na seção de testes. Salienta-se também que tanto a interface desenvolvida como também a estimativa da onda de choque a partir de 53 dados de saída do programa comercial MISSILE DATCOM® são contribuições deste desenvolvimento de tese. Os tópicos subsequentes descrevem a metodologia supracitada. 3.1 Missile Datcom Durante a década de 70 a Força Aérea dos Estados Unidos da América (USAF) reuniu todos os documentos associados a metodologias desenvolvidas para se estimar as características de estabilidade e controle de aeronaves e mísseis. Este esforço resultou em 10(dez) volumes contendo um grande número de métodos e dados experimentais utilizados por muitos anos pela comunidade aeroespacial americana. A partir desta compilação, foram estabelecidas equações, fatores de correlação e tabelas de dados experimentais suficientes para que o projetista pudesse gerar o projeto conceitual de um veículo sem grandes custos em termos de túnel de vento e/ou ensaios em voo, BLAKE (1999). Nos anos 1980 foi realizado pela McDonnell Douglas, sob contrato com o Governo dos Estados Unidos da América, um estudo de viabilidade para verificar a possibilidade de implementação de um código computacional que automatizasse os métodos compilados pela USAF para o estudo de configurações geométricas do tipo foguetes e mísseis. Tanto a estrutura dos códigos como também a forma de implementação das metodologias foram recomendados por, BLAKE (1999). O produto final deste desenvolvimento foi o código MISSILE DATCOM®. Desde então, várias versões foram disponibilizadas para o público, cada uma delas com uma contribuição de uso prático na indústria aeroespacial. A Tabela 3.1 exibe algumas das capacidades implementadas. A revisão 5/97, em particular, foi largamente difundida em muitos centros científicos por todo o mundo. Entre as particularidades desta versão destacam-se: (i) programa escrito em FORTRAN 77, (ii) além dos coeficientes de força e momento característicos de cada configuração geométrica de veículo, são disponibilizadas todas as derivadas dinâmicas necessárias para o cômputo da dinâmica do voo do veículo gerado, (iii) são disponibilizadas tanto a distribuição de pressão ao longo do corpo do veículo como também a distribuição de velocidades (em termos de 54 número de Mach), (iv) são disponibilizados em arquivos separados todas as configurações trimadas de um dado veículo, entre outros aspectos. Hoje, grande parte destas soluções foram automatizadas na forma gráfica, de forma independente por inúmeros centros de pesquisa, utilizando-se códigos como, por exemplo, C++ e MATLAB®. Este trabalho exibe uma das muitas aplicações que o código MISSILE DATCOM® viabiliza para o projetista. É interessante ressaltar que a implementação realizada neste trabalho pode ser utilizada de forma direta por desenvolvedores de códigos de engenharia, tais quais o MISSILE DATCOM® ou MISL3®. Trata-se de uma primeira aproximação para o cálculo da forma e localização de ondas de choque curvas descoladas a partir da distribuição de número de Mach que ocorre em corpos de revolução com diferentes configurações de nariz. Não é objetivo deste trabalho detalhar a metodologia matemática com que o código computacional MISSILE DATCOM® estima coeficientes aerodinâmicos de configurações geométricas do tipo míssil, foguetes ou bombas e/ou distribuições de pressão e número de Mach ao longo do corpo do veículo em estudo. Estas informações são utilizadas de forma direta, tal qual ocorre em escritórios de engenharia, e uma vez obtidos os resultados, procede-se então para a análise crítica. 3.2 Onda de choque destacada A predição da forma e localização de ondas de choque destacadas voltou a ser interesse de pesquisa em anos recentes. Este interesse tem sido estimulado pela necessidade do uso de narizes rombudos e bordos de ataque em configurações projetadas para voo hipersônico. A habilidade em predizer a forma e localização de ondas de choque destacadas é de importância primária na análise do aquecimento aerodinâmico e extremamente utilizada na escolha do tamanho máximo do modelo para um dado túnel de vento (POPE e GOIN, 1965). Toda configuração aeronáutica em regime de voo supersônico terá uma onda de choque precedendo ou colada ao seu nariz. Neste contexto, é instrutivo ressaltar a diferença encontrada na fenomenologia que aparece no campo de escoamento sobre 55 corpos que apresentam nariz ponteagudo (“sharp nose”) e rombudo (“blunt nose”). No caso de corpos rombudos, a onda curva sempre permanece destacada de forma semelhante ao mostrado na Figura 3.1. Entretanto, para qualquer corpo ponteagudo, existe um número de Mach abaixo do qual a onda choque é destacada mas acima da qual a onda choque é colada de um modo característico como uma onda Mach, conforme mostrado por BOYER (2001). Para narizes ponteagudos este número de Mach representa o limite superior da região transônica. No caso de corpos rombudos, não existe limite superior bem definido, BOYER (2001). Tabela 3. 1– Versões do código Missile Datcom Versão Revisão Melhoria acrescentada 3 12/88 Substituição de dados expandidos Incremento na configuração 4 7/89 Substituição de dados expandidos 5 4/91 Entradas em velocidades Transônicas e subsônicas, adição de arrasto Efeitos da pluma sobre o corpo Seis tipos de protuberância no corpo Modificação do centro de pressão da aleta lateral 6 6/93 Estação de trabalho UNIX, Compatibilidade com PC Flaps no bordo de fuga Aletas dobráveis Entradas semi submersas 7 5/97 Compatibilidade com Fortran 90 Expansão das derivadas dinâmicas Numerosos estudos experimentais e teóricos têm sido desenvolvidos para a determinação da forma e localização de ondas de choque bem como dos respectivos fatores de influência (KRAUSE, REINARTZ, BALLMANN, 2006). Teoricamente, a solução do campo de escoamento em torno de um corpo é dividida em duas classes distintas: aquelas em que o campo do escoamento é completamente subsônico e aquelas em que o escoamento é completamente supersônico. Para cada um destes regimes de velocidades existem métodos específicos de resolução consagrados na literatura técnica. Estes métodos são fundamentalmente dependentes da geometria em 56 estudo. Para o caso de geometrias complexas e campos de escoamentos combinados (transônicos), a literatura indica métodos mais específicos, baseados em Dinâmica dos Fluidos Computacional (CFD). De fato, há tempos, existem pesquisas associadas à formação, localização e estrutura de ondas de choque. Muitos destes estudos estão centrados em detalhes particulares importantes do problema e, por conseguinte estão restritos localmente em escopo (por exemplo, estudos restritos a velocidades hipersônicas ou a regiões próximas ao nariz). Outros estudos têm sido mais gerais e apresentam métodos para o cálculo da distância do descolamento e forma do choque sem restrições no que diz respeito a regime de velocidades. Estes métodos usualmente envolvem procedimentos iterativos laboriosos uma vez que a metodologia está baseada em CFD. Para fins de cálculo de engenharia (procedimentos rápidos de uso prático), um trabalho analítico interessante de se citar é aquele devido à MOECKEL (1949). Naquele trabalho é desenvolvido um método aproximado baseado em uma forma simplificada da aplicação da equação da continuidade para predizer a localização e forma da onda de choque à frente de corpos bidimensionais e simétricos axialmente. De forma a reduzir o escopo a um problema equivalente unidimensional, foi assumido que a forma da onda de choque entre seu ponto mais avançado e seu ponto sônico é adequadamente representada por uma hipérbole assintótica às linhas de Mach do escoamento não perturbado e a linha sônica entre a onda de choque e o corpo é representada por uma reta cuja inclinação é função do número de Mach do escoamento livre. Com estas hipóteses, a localização do choque em relação ao ponto sônico do corpo é independente da forma do nariz ou do bordo de ataque à frente do ponto sônico e se torna uma função unívoca do número de Mach. O presente trabalho de pesquisa foi realizado para obter informações rápidas sobre as características geométricas das ondas de choque destacadas. O objetivo final é desenvolver um método de engenharia baseado em soluções analíticas e resultados experimentais para prever a forma e a distância da onda de choque descolada. É importante ressaltar que este método usa dados provenientes de códigos de engenharia, tais como MISSILE DATCOM® (e MISL3®, por exemplo). Neste contexto, é possível incluir a influência da geometria do nariz do veículo na determinação das 57 características de choque destacado, tendo em vista que uma das saídas que o código computacional MISSILE DATCOM® disponibiliza é a distribuição de número de Mach na superfície do corpo em análise. 3.3 Formulação Matemática Este tópico apresenta a formulação matemática adotada para estimar a distância e a forma da onda de choque que se forma em um corpo que se move a velocidades supersônicas a zero grau de ângulo de ataque. O campo de escoamento a ser considerado é mostrado na Figura 3.1. A metodologia se baseia no trabalho de MOECKEL (1949), com a diferença de que nesta pesquisa é considerada a geometria do nariz no cômputo da forma e localização da onda de choque. Essencialmente, a formulação matemática é baseada em três modelos, quais sejam: geometria da onda de choque, princípio da continuidade e distância entre a onda de choque e o corpo. Figura 3. 1– Geometria do problema, MOECKEL (1949) As hipóteses adotadas na resolução do problema estão enumeradas abaixo: 58 i. O choque curvo é representado por uma hipérbole; ii. O ângulo �, (Figura 3.1), é uma função da deflexão da linha de corrente no ponto sônico da onda de choque e da deflexão da linha de corrente no ponto sônico sobre o corpo; iii. É adotada a analogia entre o escoamento unidimensional interno e o escoamento com ondas de choque destacadas. A linha sônica é considerada normal na vizinhança da direção do escoamento médio; iv. A influência do nariz do corpo na forma da onda de choque e na distância entre o corpo e a onda de choque é obtida a partir da distribuição de Mach (usando software comercial) ao longo da superfície do corpo. 3.4 Geometria da onda de choque Algumas características são típicas das ondas de choque destacadas (BOYER, 2001; LOVE, 1957 e MOECKEL, 1949): i. São normais ao escoamento não perturbado no ponto mais próximo do nariz do corpo; e ii. São assintóticas ao escoamento não perturbado a grandes distâncias do ponto mais próximo ao nariz do corpo. Uma curva simples que possui estas características é a hipérbole representada pela equação (3.1): � 2 0 2 xx y � � , (3.1) sendo � � � ! " �� � � �� � � � � 0 1 1sincot M � , 0M é o número de Mach do escoamento livre e 0x é a distância entre o vértice da onda e a interseção de suas assíntotas. Com 59 esta hipótese, o ângulo entre a direção do escoamento e a tangente ao choque em qualquer ponto é dada pela Eq.(3.2): y xy dx dy 2 2 0 22 tan � � . � �� , (3.2) A posição ( Sx , Sy ) do ponto sônico na curva da onda de choque é dada por: S S x x .� � 22 0 cot� � (3.3) S S S x y .�� . 22 0 cot cot � � (3.4) sendo S. é a inclinação da curva no ponto sônico. Se for usada a coordenada y do ponto sônico do corpo ( SBy ) como dimensão de referência, então a forma da onda de choque pode ser modelada conforme a Eq. (3.5): 2 0 2 1 �� � � �� � � ��� � � �� � � � SBSBSB y x y x y y � (3.5) Da definição de Sy é possível concluir que: 1tan220 �� S SB S SB y y y x .�� . (3.6) Neste contexto, a localização adimensional do ponto sônico na onda de choque é dada por: S SB SB S y x y x .� � 22 0 cot� � . (3.7) 60 3.5 Distância entre o choque e o nariz do corpo A partir da Figura 3.1 é possível concluir que a distância entre o ponto mais avançado do choque ( 0x ) e a coordenada x do ponto sônico ( SBx ) está relacionada com a distância entre o choque e o corpo (MOECKEL, 1949) de acordo com a equação (3.8) : SBSB SB SB y x y x y L 0�� . (3.8) Mas �tan1�� � � �� � � ��� SB S SB S SB SB y y y x y x , (3.9) então: � � �� tantan ��� C y y y L SB S SB , (3.10) sendo � �1tantan 22 ��� SSC .�.�� (3.11) e 2 SBS �� � � � (3.12) Note que a variável S� representa a deflexão da linha de corrente no ponto sônico da onda de choque e SB� é a deflexão da linha de corrente no ponto sônico do corpo. A expressão apresentada na Eq. (3.12) é a principal hipótese utilizada no trabalho de MOECKEL (1949). 61 3.6 Princípio da continuidade A fim de se estimar a quantidade SB S y y , a relação da continuidade é aplicada para o fluido que escoa em torno do corpo. A forma integral da equação de continuidade é inútil para o presente método, desde que, para o seu uso, é necessário conhecer, a priori, a distribuição de pressão e velocidade do escoamento ao longo da linha sônica. O princípio da continuidade é definido pela analogia entre o escoamento unidimensional interno e o escoamento com ondas de choque destacadas, a linha sônica é considerada normal na vizinhança da direção do escoamento médio. A partir desta hipótese: � � � �1,,,, 000000 ���� SSSSS MTTCdAPMTCdAPm ��� , (3.13) onde P é a pressão de estagnação, T é a temperatura de estagnação, A é a área de seção transversal e 0Cd é um coeficiente apresentado na Eq. (2.9). O subscrito “0” representa a condição de escoamento livre e o subscrito “S” está relacionado com o ponto sônico na onda de choque. Desde que a pressão de estagnação permanece constante ao longo de cada linha de corrente antes do choque, adota-se a hipótese que estabelece que há um valor médio para a pressão de estagnação ( Average SP ) o qual permanecerá inalterado entre a onda de choque e a linha sônica. Desde que a temperatura total é constante, a equação da continuidade simplificada pode ser escrita como: � � � �12 1 2 00 0 0 1 12 � � � � � � ! " � �� � � � � � M P P M A A Average S S . (3.14) Usando a relação de área na Figura 3.1é possível escrever: 62 � � 2 1 cos1 ��� �B y y SB S , (3.15) sendo: � � � �12 1 2 00 0 1 12 � � � � � � ! " � �� � � � � � M P P MB Average S (3.16) Neste contexto, a distância do choque é dada por: � � � � ��� tantancos1 2 1 ���� � CB y L SB . (3.17) 3.7 Algoritmo Numérico Essencialmente, a onda de choque divide o campo em duas partes. Em toda parte a montante da onda de choque, o escoamento é uniforme e a pressão total ou de estagnação é constante. A jusante da onda de choque, o escoamento varia em todo o campo e cada linha de corrente tem uma pressão de estagnação diferente. Esta variação da pressão de estagnação é devido à variação na entropia através da onda de choque. A variação de entropia é função do número Mach do escoamento livre e do ângulo de desvio ( S� ) da linha de corrente. A metodologia desenvolvida por MOECKEL (1949), usa um valor médio adequado para a pressão de estagnação no cômputo da vazão mássica em torno do corpo. Assume- se que esta pressão de estagnação média é calculada como sendo aquela associada ao centroide de massa do fluido localizado entre o ponto de estagnação e a linha sônica. A linha de corrente que passa pelo centroide entra na onda de choque em 2 S C y y � para escoamento plano e em 3 2 S C y y � para escoamento axialmente simétrico. Neste trabalho é desenvolvida uma metodologia que não usa esta aproximação. É usado um algoritmo numérico com o qual é possível 63 estimar a forma e a localização da onda de choque através de um processo iterativo. Os parágrafos seguintes descrevem o algoritmo numérico: i. Dados de entrada: Condição do escoamento livre ( 000 ,, TPM ), geometria do corpo: � �SBSBSB xyy � e distribuição do número de Mach sobre o corpo: � �SBSBSBSB yxMM ,� ; ii. Determinar � �1�� SBSBSB Myy : localização sobre o corpo, onde o número de Mach vale 1; iii. Determinação do ângulo ( S. ) da onda de choque no ponto sônico: 1 sin 2 11 cos 2 1sin 2 11 22 0 22 0 22 0 2 0 2 � � � � � � � � � S S S After M M M M M .� . �.� � . (3.18) iv. Determinar � . 2 SBSk �� � � � (3.19) A priori, é assumido o valor de 1 (um) para o parâmetro k. Este parâmetro deve ser identificado durante o processo de comparação com os dados experimentais. Note que S� pode ser obtido a partir da equação da onda de choque oblíqua, qual seja: � �1sin 2 1 1sincottan 22 0 2 0 22 0 �� � � � S S SS MM M .� ..� . (3.20) v. Determinar 0x , Eq. (3.6). vi. Determinar as coordenadas do ponto sônico na onda de choque (Eq. 3.3 e 3.4) vii. Determinar L, Eq.(3.10) 64 viii. Determinar SBx . Da Figura 3.1: LxxSB �� 0 . (3.21) ix. Determinar New� da geometria do problema: SBS SSBNew yy xx � � �� . (3.22) x. Verificar a convergência do método: � �� / � � New . (3.23) Se � � 310�0/Abs , vá para o item iv. 3.8 Método de Identificação de Sistemas O algoritmo apresentado anteriormente, leva em consideração a hipótese de MOECKEL (1949), qual seja 2 SBS �� � � � . Ressalta-se que este parâmetro foi obtido a partir da observação de resultados experimentais gerados para algumas configurações de nariz. A partir das simulações realizadas por este trabalho foi possível verificar que tanto a forma da onda de choque como também a distância a partir do corpo é bastante sensível à variação de �. Neste contexto, procedeu-se no sentido de se desenvolver uma metodologia que atualize o valor do referido ângulo a partir de resultados experimentais. A técnica empregada é também denominada Identificação de Parâmetros (“Systematic Identification Parametric – SIP”). Essencialmente, esta técnica consiste em incluir um otimizador que poderá ajustar os coeficientes desconhecidos(neste caso o ângulo �) ou variáveis até se obter uma aproximação melhor do ajuste aos dados experimentais usando a técnica dos mínimos quadrados. O SIP tem várias vantagens sobre outra opção de modelagem como, por exemplo CFD. Criar um modelo matemático que captura as características de ondas de choque usando CFD é viável. No entanto, é 65 extremamente caro em termos de custo computacional. O objetivo deste trabalho é desenvolver um método de engenharia que proporcione resultados rápidos e de baixo custo computacional com precisão azoável. O SIP ignora a complexidade do campo de escoamento (ondas de choque, camada limite e assim por diante) e foca sobre os resultados finais. Isso faz com que o método SIP não seja apenas mais simples, mas também rápido o suficiente para produzir resultados de engenharia. Outro aspecto interessante desta metodologia é que o SIP também lida com condições não lineares, como, por exemplo, transferência de calor, presença de descontinuidades no campo de escoamento, entre outros, viabilizando o seu emprego em uma extensa faixa de problemas de engenharia. O procedimento SIP adotado neste trabalho segue a seguinte linha de raciocínio: em primeiro lugar, dados físicos e geométricos são lidos de um arquivo de dados. A partir destas informações é possível determinar a configuração da onda de choque (forma e distância do corpo) usando o algoritmo numérico acima mencionado. O ângulo descrito pela Eq.(3.19) é uma primeira abordagem, e, portanto, é necessário atualizá-la. Neste contexto, a variável k é enviada para o otimizador. Neste trabalho é utilizado o otimizador fmincon, o qual é disponibilizado pelo código computacional MATLAB®. Este otimizador utiliza o método de Levenberg-Marquardt sujeito a restrições definidas pelo usuário. A função do otimizador chama uma função implementada pelo usuário a qual calcula a forma e a localização da onda de choque com os parâmetros atuais fornecidos pelo otimizador. Após o cálculo das características principais da onda de choque, é retornada uma função de erro para o otimizador, baseada nas diferenças entre a distância de choque previsto e experimental. Na sequencia, o otimizador ajusta os parâmetros desconhecidos (neste caso �, através do parâmetro k) e o processo se repete até que determinados critérios de convergência sejam atendidos. É importante salientar que o otimizador lsqcurvefit não funcionou muito bem em escoamento com alto número de Mach. A estimativa para o ângulo � não foi realista, quando comparada com a literatura de aerodinâmica. 66 4. DIMENSIONAMENTO DE UM TVS 4.1 Requisitos O anteprojeto de um túnel TVS exige um trabalho multidisciplinar bastante elaborado. Neste trabalho estão envolvidos o aerodinamicista, o engenheiro de instrumentação e controle, o especialista em acústica, o engenheiro de projeto mecânico e o engenheiro de estrutura. Apesar da relativa simplicidade de construção de TVS, os problemas associados à relação entre a qualidade de escoamento na seção de testes e o sistema de controle desafiam a competência e a habilidade de seus projetistas (PAGLIONE, 1978). Neste tópico do trabalho são elencados e discutidos os principais requisitos que devem ser adotados na concepção de um TVS. Estes requisitos são fundamentalmente dependentes do tipo de ensaio (modelo, tempo de duração dos ensaios, número de Reynolds e número de Mach na seção de testes) que se pretende realizar. A Tabela 4.1 exibe os requisitos básicos utilizados no projeto do modelo para túnel de vento. A Tabela 4.2 indica os requisitos básicos utilizados na concepção geométrica e operacional de um TVS. Ressalta-se que a razão máxima de bloqueio do modelo é definida como sendo a relação entre a maior área transversal do modelo e a área da seção de ensaio. Os valores numéricos indicados ilustram os requisitos adotados por um túnel de vento TVS acadêmico. Tabela 4. 1 – Requisitos para concepção do modelo REQUISITO VALOR Razão de bloqueio # 1 % Dimensões máximas do modelo Corda: 10 cm Largura: 15 cm 67 Tabela 4. 2 – Requisitos para concepção do TVS REQUISITO VALOR Faixa de número de Mach ( STM ) 2,0 – 4,0 Tempo de corrida (mínimo) ( runt ) 60 s (Mach 2,0 – 3,5) 40 s (Mach 4,0) Número de Reynolds ( yRe ) Comprimento de referência: lado da seção de testes Coeficiente de expansão politrópica (n) 2,11n (Processo politrópico no reservatório de ar comprimido) Válvula reguladora de pressão Curva de vazão mássica Perda de pressão entre reservatório de ar e câmara de Tranquilização Nível de Turbulência na Câmara de Tranquilização 0,5 % a 1,0 % Velocidade média na Câmara de