Embed Size (px)
Citation preview
Comportamento Dinâmico de Corpos no Seio de um Fluido
Vitor Fernando Rosa Caetano
Dissertação para obtenção do Grau de Mestre em
Engenharia Mecânica
Júri
Presidente: Prof. Doutor Luís Rego da Cunha de Eça
Orientador: Prof. Doutor José Carlos Fernandes Pereira
Co-Orientador: Doutor José Manuel da Silva Chaves Ribeiro Pereira
Vogal: Prof. Doutor João Manuel Gonçalves de Sousa Oliveira
Dezembro 2010
ii
Agradecimentos
Ao terminar esta Tese de mestrado registo os meus sinceros agradecimentos às pessoas
que de várias formas contribuíram para que esta fosse uma realidade.
Em primeiro lugar quero agradecer Professor José Carlos Pereira e ao Doutor José Chaves
Pereira, por me possibilitarem realizar esta dissertação, bem como pela disponibilidade e amizade
que sempre mostraram.
Aos meus pais e à minha irmã, quero aqui expressar o meu profundo agradecimento por
todo o apoio que sempre me deram durante toda a minha vida.
Agradeço ao amigo Filipe Garrido, pelo apoio e entusiasmo que sempre demonstrou ao longo
do período em que este trabalho se realizou.
Por fim, mas não com menor importância, agradeço à minha namorada, Margarida
Rodrigues, pelo conforto, compreensão e apoio que incondicionalmente sempre demonstrou.
iii
Resumo
Nesta Tese, aplicou-se a mecânica dos fluidos computacional à dinâmica de corpos que se
movem no seio de um fluido. Foi desenvolvida uma rotina em FORTRAN para modelar a dinâmica do
corpo rígido com seis graus de liberdade, recorrendo a referenciais não inerciais. Foi feito o
acoplamento dessa rotina a um programa comercial (Star-CD), de cálculo de escoamentos de
fluidos. Foram também trabalhados diversos assuntos como a interacção desse programa com
rotinas externas e a técnica das malhas móveis, incluindo processos de computação em paralelo.
Seguidamente foram realizadas uma série de simulações bidimensionais e tridimensionais
por forma a dar validade ao modelo acoplado. Por fim utilizou-se esse modelo para chegar a alguns
resultados relevantes sobre a queda livre de placas leves na atmosfera (falling plates), chegando-se
a estudar algumas tipologias de movimentos típicos conhecidos como o tumbling.
Palavras Chave: Mecânica dos Fluidos Computacional; Dinâmica do Corpo Rígido;
Referenciais não Inerciais; Falling Plates; Tumbling; Fluttering.
iv
Abstract
In this Thesis Computational Fluid Dynamics (CFD) was applied together with the
Computational Rigid Body Dynamics (RBD) to the study of moving bodies immersed in a fluid
environment. A Fortran subroutine using non-inertial reference frames was developed to model the 6
degrees-of-freedom (6dof) rigid body dynamics. This routine was coupled to the commercial CFD
software package Star-CD, using the available Star-CD subroutine interfaces. All the interactions
between the two physical models as well as the moving mesh techniques were fully tested and
validated, including with parallel processing. The validation was successful for both 2D and 3D
models.
A set of simulations was conducted to extensively analyze the 2D motion of falling plates
under several different conditions. The tumbling and fluttering movement patterns were identified to
appear in different flow regimes, in accordance with other results found in scientific literature.
The employed strategies led to a process suitable for general 3D simulations of moving rigid
bodies interacting with a fluid environment.
Keywords: Computational Fluid Dynamics; Rigid Body Dynamics; Non-inertial reference frames;
Falling Plates; Tumbling; Fluttering.
v
Índice
Agradecimentos ...................................................................................................................................... ii
Resumo ...................................................................................................................................................iii
Palavras Chave: ..................................................................................................................................iii
Abstract .................................................................................................................................................. iv
Keywords:........................................................................................................................................... iv
Nomenclatura ......................................................................................................................................... ix
1. Introdução ........................................................................................................................................ 1
1.1. Motivação ................................................................................................................................ 1
1.2. Objectivo ................................................................................................................................. 1
1.3. Revisão Bibliográfica. .............................................................................................................. 2
1.4. Interacção Fluido-Estrutura ..................................................................................................... 3
1.5. Tumbling e Fluttering .............................................................................................................. 4
2. Modelação da Dinâmica do Corpo Rígido ....................................................................................... 5
2.1. Sistemas de Coordenadas ...................................................................................................... 5
2.2. Derivada em Ordem ao Tempo de um Vector Escrito num Referencial Animado de
Movimento de Rotação ....................................................................................................................... 6
2.3. Dinâmica Tridimensional do Corpo Rígido ............................................................................. 8
2.3.1. Dinâmica da Translação ................................................................................................. 8
2.3.2. Dinâmica da Rotação - Equações de Euler .................................................................... 9
2.4. Matriz de Rotação e Ângulos de Euler. ................................................................................ 12
2.5. Rotação do Corpo Rígido – Derivada da Matriz Rotação ..................................................... 14
3. Modelação do Escoamento ........................................................................................................... 16
3.1. Equações Base ..................................................................................................................... 16
3.2 Método de Discretização ................................................................................................... 17
3.2.1 Esquemas de Discretização do Termo Convectivo .......................................................... 17
3.2.2 Método de Acoplamento Pressão Velocidade .................................................................. 19
3.2.3 Esquemas de Discretização do Termo Temporal ............................................................. 19
3.3. Geração de Malha ................................................................................................................. 21
3.3.1. Condições de Fronteira Utilizadas ................................................................................ 21
vi
3.3.2 Mecanismos de Geração de Malha .................................................................................. 22
3.3.3. Malha 2D ....................................................................................................................... 24
3.3.4. Malha 3D ....................................................................................................................... 30
4. Algoritmo de Cálculo ...................................................................................................................... 33
5. Resultados ..................................................................................................................................... 36
5.1. Simulações Preliminares ....................................................................................................... 36
5.2. Resultados das Modelações Bidimensionais (2D)................................................................ 37
5.2.1 Parâmetros Adimensionais Utilizados na Apresentação dos Resultados ........................ 38
5.2.2. Validação ....................................................................................................................... 41
5.2.3. Relação Entre Velocidades e Número de Reynolds. .................................................... 45
5.2.4. Simulação com Movimento de fluttering ...................................................................... 48
5.2.5. Estudo das Trajectórias em Função da Razão de Densidades .................................... 49
5.2.6. Estudo das Trajectórias em Função do Ângulo de Lançamento Inicial ........................ 51
5.2.7. Campos de velocidade, de pressão estática e vórticidade. .......................................... 55
5.2.8. Outros Resultados Ilustrativos de Simulações 2D ........................................................ 59
5.3. Resultados de Uma Simulação Tridimensional (3D) ............................................................ 61
6. Conclusões .................................................................................................................................... 71
Referências ........................................................................................................................................... 72
vii
Índice de Figuras
Figura 1 - As várias topologias de trajectórias de placas em queda no seio de um fluido (Fonte
[11]). 4
Figura 2 - Sistemas de coordenadas 6
Figura 3 - Ângulos de Euler 12
Figura 4 - Molécula Computacional para o Esquema Implícito (Fonte: [21]) 20
Figura 5 - Método de Geração de Malha Automática 23
Figura 6 - A Malha; B Cotagem da Malha 24
Figura 7 - Pormenores da Parte Central da Malha 2D 25
Figura 8 - Método de Elaboração da Malha 2D 25
Figura 9 - Pormenor das Células Junto à Placa 26
Figura 10 - Pormenor da Parte Fixa da Malha 27
Figura 11 - Pormenor da União da Malha Móvel e Malha Fixa 28
Figura 12 – Várias Posições da Parte da Malha com Capacidade de Rotação 28
Figura 13 - Condições Fronteira da Malha 2D: A - Attach; B - Entrada, Saída e Semi-Plano; C -
Parede 29
Figura 14 - Domínio da Malha Tridimensional 30
Figura 15 - Cortes da Malha 3D 30
Figura 16 - Parte Interior da Malha 31
Figura 17 - Pormenor da Malha em Volta do Disco 32
Figura 18 - Algoritmo de Cálculo 33
Figura 19 - Rotinas Externas 35
Figura 20 – Comparação de Trajectórias 41
Figura 21 – Força na Direcção y 42
Figura 22 – Força na Direcção x 43
Figura 23 - Velocidades 44
Figura 24 – Velocidade Angular 44
Figura 25 - Resultados para e 45
Figura 26 – Resultados para e 46
Figura 27 - Resultados para α=8 e =1198 47
Figura 28 - Trajectória do Movimento Fluttering 48
Figura 29 - Componentes das Forças Adimensionalizadas 49
Figura 30 - Velocidade Angular Adimensional 49
Figura 31 - Trajectórias em Função de B 50
Figura 32 - Pormenor Inicial das Trajectórias em Função de B 51
Figura 33 - Trajectórias em Função do Ângulo Inicial 52
Figura 34 - Pormenor das Trajectórias em Função do Ângulo Inicial 53
Figura 35 - Trajectórias com Velocidade Inicial Nula 54
viii
Figura 36 - Campos de Velocidades Adimensionais. Da esquerda para a direita e de cima para
baixo, os instantes (0.44; 0.47; 0.50; 0.53; 0.56; 0.59; 0.62; 0.65) s 56
Figura 37- Campos de Pressões Adimensionais. Da esquerda para a direita e de cima para
baixo, os instantes (0.44; 0.47; 0.50; 0.53; 0.56; 0.59; 0.62; 0.65) s 57
Figura 38 - Campos de Vorticidade. Da esquerda para a direita e de cima para baixo, os
instantes, (0.44;
0.47; 0.50; 0.53; 0.56; 0.59; 0.62; 0.65) s 58
Figura 39 - Pormenor de um Campo de Velocidades 59
Figura 40 - Vectores Velocidade Junto à Placa 59
Figura 41 - Alguns Exemplos de Trajectórias. Da esquerda para a direita e de cima para baixo
a, b, c e d 60
Figura 42 - Referencial Utilizado. Posição Inicial do Disco. 62
Figura 43 - Trajectória 3D. Nota: estes eixos não estão normalizados 62
Figura 44 - Trajectórias nos Planos xoy e yoz. Nota: Chama-se a atenção para o facto de a
escala não ser a mesma nos dois eixos. 63
Figura 45 - Trajectórias nos Planos xoz. Nota: Chama-se a atenção para o facto de a escala
não ser a mesma nos dois eixos. 63
Figura 46 - Componentes da Forças Instantânea que Actua no Disco 64
Figura 47 - Componentes dos Momentos Instantâneos que Actuam no Disco 65
Figura 48 - Componentes da Velocidade do Disco 65
Figura 49 - Velocidade Angular do Disco 66
Figura 50 - Ângulos de Euler do Disco 66
Figura 51 - Campos de Velocidades em m/s. Da esquerda para a direita e de cima para baixo,
os instantes (0.71; 0.74; 0.77; 0.80; 0.83; 0.86; 0.89; 0.92) s 67
Figura 52 - Campos de Pressão em Pa. Da esquerda para a direita e de cima para baixo, os
instantes (0.71; 0.74; 0.77; 0.80; 0.83; 0.86; 0.89; 0.92) s 68
Figura 53 - Campos de Vorticidade. Da esquerda para a direita e de cima para baixo, os
instantes (0.71; 0.74; 0.77; 0.80; 0.83; 0.86; 0.89; 0.92) s 69
Figura 54 - Pormenor de um Campo de Velocidades 70
Figura 55 - Vectores Velocidade 70
ix
Nomenclatura
Aceleração
Distância
, , Vectores unitários nas direcções dos eixos coordenados, versores
Força
Aceleração da gravidade
Espessura
G Origem do sistema de coordenadas; centro de gravidade
Gxyz Sistema de coordenadas rectangulares com translação e rotação
GX'Y'Z' Sistema de coordenadas rectangulares com orientação fixa
Momento de inércia
Comprimento
Momento angular
m Massa
Momento de força
O Origem das coordenadas
Oxyz Sistema de coordenadas rectangulares com translação e rotação
OXYZ Sistema de coordenadas rectangulares absoluto
Pressão
Matriz de rotação
t Tempo
Intervalo de tempo
Campo de velocidades em cada instante
Velocidade
Vector posição
x,y,z Coordenadas rectangulares; distâncias
X,Y,Z Coordenadas rectangulares
X'Y'Z' Coordenadas rectangulares
Velocidade angular
Viscosidade cinemática do fluido
Viscosidade dinâmica
Densidade
Ângulos de Euler
1
1. Introdução
1.1. Motivação
O comportamento dinâmico de corpos no seio de fluidos é um tema de capital relevância
para várias aplicações, pois esta área do conhecimento descortina os fenómenos físicos que estão
na base de alguns problemas práticos, tais como: a largada de bombas ou corpos para alijar uma
aeronave; a projecção na atmosfera de corpos vários, como pequenos troncos ou ramos
eventualmente em combustão no caso dos incêndios florestais [1]; lançamento de projécteis de
artilharia [2]; ou corpos de utilidade lúdica ou desportiva; passando pela dinâmica de partículas de
reduzidas dimensões.
O assunto que se pretende estudar também está relacionado com a temática do voo de
insectos [3,4,5]. A evolução do conhecimento nesta área também pode ter interesse em fenómenos
como certos tipos de sedimentação, dispersão de sementes e em algumas técnicas de separação de
materiais em engenharia química.
1.2. Objectivo
O objectivo deste trabalho inclui-se no domínio da mecânica dos fluidos computacional
aplicada à dinâmica de corpos que se movem no seio de um fluido sujeitos à força da gravidade.
Pretende-se obter uma modelação matemática para o fenómeno, na situação bidimensional e
também tridimensional, que será realizada recorrendo ao código comercial Star-CD, de grande
versatilidade a nível de geração de malha e variedade de opções físicas e numéricas. A modelação
do movimento do corpo será realizada com recurso a referenciais não inerciais e à dinâmica do corpo
rígido com seis graus de liberdade, a qual será acoplada à dinâmica de fluidos através do uso de
rotinas programadas em FORTRAN especialmente desenvolvidas para este tipo de problemas. É
também objectivo deste trabalho obter resultados credíveis, susceptíveis de serem comparados com
os resultados de outros trabalhos sobre o assunto, numéricos e experimentais. Tais resultados
incluem trajectórias do centro de massa, forças instantâneas, velocidade linear e velocidade angular
instantânea.
2
1.3. Revisão Bibliográfica
Folhas, sementes de árvores e cartões de papel que oscilam para o lado e pairam (fluttering),
ou rodam sobre si próprios e deslizam para um dos lados (tumbling), são maravilhosos exemplos do
quotidiano, de objectos sólidos que se movem no seio de um fluido. Para prever as suas complexas
trajectórias é necessário o conhecimento das forças instantâneas que actuam no sólido. Devido à
complexidade do fenómeno, existem muito poucos resultados analíticos sobre o assunto, e os que
existem são baseados em simplificações.
A dança graciosa da queda de uma folha há muito que inspira tanto físicos como poetas, pois
o estudo de placas em queda livre tem já uma longa história, já em1854, o teórico James Clerk
Maxwell descreveu qualitativamente o fenómeno denominado por tumbling [6], ainda antes do
desenvolvimento da aerodinâmica clássica. Desde então alguns cientistas esporadicamente têm
tentado explicar o fenómeno dos objectos planos a cair na atmosfera. Entre 1950 e 1960 houve uma
explosão de interesse sobre o assunto principalmente por motivos militares. No entanto nas últimas
décadas, a maioria das investigações apenas apresentam resultados qualitativos ou em termos de
propriedades médias. Como foi o caso de A. Belmonte, H. Eisenberg and E. Moses, em 1998, que
estudaram a transição de fluttering para tumbling, variando o momento de inércia adimensionalizado
em escoamentos experimentais quasi-bidimensionais [7].
Ao analisar o movimento da queda de um cartão de papel, Z. Jane Wang e seus
colaboradores, Anders Anderson e Umberto Pesavento, em 2004 e 2005, usando vídeo digital de alta
velocidade conseguiu medir com precisão a trajectória da queda do papel e também da queda de
placas de alumínio em tanques com água. Este método tem resolução suficiente para determinar de
imediato as acelerações, a partir das quais se podem calcular as forças e momentos instantâneos do
movimento da queda das placas. Tais experiências com chapas de alumínio em água foram
concebidas de forma a minimizarem os efeitos tridimensionais. Posteriormente, usando as
potencialidades da mecânica dos fluidos computacional, A. Andersen, U. Pesavento e Z. Jane
Wang, obtiveram soluções numéricas da equação de Navier-Stokes para escoamento bidimensional
tendo por base uma teoria invíscida. [8,9,10].
David L. Finn (2007) estudou as trajectórias de cartões em queda livre, como um modelo
dependente da razão de densidades entre o corpo e o fluido e da razão entre o comprimento e a
espessura do corpo [11].
Mais recentemente, em 2008, Changqiu Jin e Kun Xu, realizaram um estudo numérico
utilizando teoria cinética de gases, tanto para placas elípticas como para rectangulares [12].
3
1.4. Interacção Fluido-Estrutura
Estamos perante um fenómeno denominado de interacção fluido-estrutura, sempre que se
pretende estudar um sistema em que um corpo rígido, ou mais ou menos rígido, interage com um
meio contínuo. Podem-se dividir estes fenómenos em três tipos essenciais:
Situações em que o corpo sólido tem movimento imposto, e que não depende
inteiramente da acção do fluido. Como por exemplo o pistão de um motor de
combustão interna.
Casos em que o corpo tem alguma flexibilidade e é actuado pela tensão do fluido,
fenómeno conhecido como aero-elasticidade, como são exemplos problemas que
envolvem membranas, turbinas eólicas, estudo de asas de micro-veículos [3,4,5,13].
Problemas em que geralmente se considera que o corpo é perfeitamente rígido, mas
o movimento deste é imposto pelo fluido e pela força da gravidade [9]. É neste grupo
que este trabalho se insere, uma vez que nestas simulações não existe qualquer
movimento imposto aos corpos. O facto de o fluido impor o movimento do corpo e
vice versa, dá a estes problemas uma dificuldade adicional, pois em geral estes
casos são difíceis de resolver, principalmente se a razão entre as massas especificas
do corpo e do fluido for próxima da unidade.
4
1.5. Tumbling e Fluttering
Quando um cartão cai livre na atmosfera, geralmente descreve movimentos graciosos. Se
este oscilar de um lado para outro o movimento denomina-se de fluttering (Figura 1-a), se por sua
vez rodar sobre si próprio e deslizar para o lado, voltando a rodar e a deslizar para o lado e assim
sucessivamente, muitas vezes invertendo o sentido descendente levantando-se contra a gravidade,
designa-se de tumbling (Figura 1-b). No entanto, existem situações em que o movimento descrito
pelo corpo é ainda mais complexo, como é o exemplo do tumbling de duas oscilações por período
(Figura 1-c). Também pode acontecer que o corpo adquira um movimento constituido por uma
mistura bem sincronizada de tumbling e fluttering, (caso da Figura 1-d). Alguns, especialmente
quando a razão entre as massas especificas do corpo e do fluido são maiores, rodam determinado
ângulo alternadamente para um lado e para outro, descrevendo trajectórias quase verticais tendo um
movimento com características do fluttering (Figura 1-f). Existem ainda outras variantes entre eles,
alguns com aparência caótica (Figura 1-e) e ainda outros que têm aparência caótica até determinado
momento, entrando depois num regime perfeitamente organizado [11,14,15].
Figura 1 - As várias topologias de trajectórias de placas em queda no seio de um fluido (Fonte [11]).
5
2. Modelação da Dinâmica do Corpo
Rígido
2.1. Sistemas de Coordenadas
Em Mecânica, um corpo rígido, também chamado corpo sólido e por vezes apenas objecto,
pode-se definir como um sistema de pontos materiais cuja distância mútua é invariável. De facto, os
sistemas que existem na natureza só podem satisfazer esta condição de forma aproximada. No
entanto, a maioria dos corpos sólidos nas condições habituais variam a sua forma e dimensões tão
pouco, que no estudo das leis do movimento do corpo rígido, considerado como um todo, podem-se
desprezar completamente estas variações.
Para descrever o movimento de um corpo rígido, consideraram-se três sistemas de coordenadas:
Um sistema de coordenadas considerado imóvel, o referencial OXYZ, o único dos três
referenciais que é inercial, com o objectivo de descrever o movimento de translação do
sólido;
O referencial GX'Y'Z', não inercial, cujos eixos se mantêm sempre paralelos com os do
referencial OXYZ, este referencial mantém-se sempre no centro de gravidade do corpo e é
em relação a este que o corpo rígido tem movimento de rotação, portanto é um referencial
de orientação fixa.
Por último, o referencial não inercial Gxyz, que se pressupõe estar vinculado rigidamente
com o corpo e participa em todos os seus movimentos. A origem deste referencial estará
sempre no centro de inércia do corpo.
6
Figura 2 - Sistemas de coordenadas
2.2. Derivada em Ordem ao Tempo de um Vector Escrito num
Referencial Animado de Movimento de Rotação
Considere-se dois referenciais, um fixo OXYZ e outro móvel GX'Y'Z', mas com orientação fixa
em relação ao primeiro e animado com uma velocidade de translação arbitraria . Suponha-se
também que um ponto material P com movimento arbitrário se desloca em relação a estes
referenciais com velocidade em relação ao referencial móvel.
(2.1)
(2.2)
Como se sabe, as componentes de um dado vector velocidade dependem do referencial
que se escolhe para o representar, mas as suas derivadas temporais obtêm-se simplesmente
7
derivando as componentes escalares para cada direcção, uma vez que os vectores unitários das
direcções dos eixos coordenados são invariantes no tempo, pressupondo um referencial de
orientação fixa. Utilizando notação indicial, com , o vector velocidade e a sua derivada
temporal tomam a forma seguinte.
(2.3)
(2.4)
Tal simplicidade não sucede se um referencial móvel estiver animado de movimento de
rotação. Para se analisar esta situação considere-se um outro referencial Gxyz com movimento de
rotação em relação a GX'Y'Z' com velocidade angular , mas com a mesma origem (Figura 3). Neste
caso tanto o módulo como a direcção de variam. Em notação indicial:
(2.5)
As derivadas em ordem ao tempo dos vectores unitários (versores) das direcções justificam-
se porque a sua direcção varia ao longo do tempo, embora a sua intensidade se mantenha
constante, uma vez que se considerou que este referencial tem movimento de rotação. Tomando
como exemplo o vector unitário e com o objectivo de calcular a sua derivada temporal, designe-se
por um ponto colocado na sua extremidade, deste modo a sua derivada temporal coincide com a
velocidade desse ponto .
(2.6)
De forma genérica para os outros eixos e em notação indicial vem,
(2.7)
(2.8)
Representando simbolicamente por
aparece o seguinte resultado, que é válido no
caso em que representa velocidade ou outra função vectorial dependente do tempo.
8
(2.9)
em que
representa a variação do vector por unidade de tempo supondo que o referencial não roda
(variação das suas componentes) e representa a contribuição da rotação do referencial para a
variação do vector. Para o caso de se pretender uma análise mais aprofundada sobre o assunto,
este pode ser encontrado mais detalhadamente em alguns manuais de Mecânica Clássica como
exemplo deixa-se a referência [16].
2.3. Dinâmica Tridimensional do Corpo Rígido
2.3.1. Dinâmica da Translação
Nesta Tese o movimento de translação realizado pelos sólidos é representado com três
graus de liberdade, que são dados pelo vector posição do centro do referencial solidário com o corpo
, em relação ao referencial inercial fixo . A modelação do movimento de translação é
realizada recorrendo à segunda lei de Newton, que nos permite saber a aceleração do corpo partindo
da força resultante que o actua:
(2.10)
onde representa o vector resultante das forças exteriores aplicado ao corpo, representa a massa
do corpo e é o vector aceleração. Sabendo que a aceleração e a velocidade do centro de
gravidade do corpo são dadas por:
(2.11)
(2.12)
(2.13)
9
pode-se determinar a sua velocidade e posição. Tal processo foi numericamente realizado
procedendo-se à discretização da derivada temporal com o método de Euler, neste caso explícito de
primeira ordem, nos instantes e ( , para um intervalo de tempo constante tal que:
(2.14)
assim, obtêm-se estimativas numéricas para a posição e velocidade do corpo no instante ( .
(2.15)
(2.16)
2.3.2. Dinâmica da Rotação - Equações de Euler
No que diz respeito à modelação da dinâmica da rotação, as equações do movimento de um
corpo rígido em rotação em torno de um ponto fixo podem obter-se recorrendo ao teorema do
momento angular, ou seja a segunda lei de Newton para a rotação:
(2.17)
onde representa o momento angular do corpo em relação ao seu centro de gravidade e
representa a resultante da soma dos momentos exteriores ao corpo. Por outro lado o momento
angular relaciona-se com a velocidade angular do sólido através da relação:
(2.18)
(2.19)
10
onde representa o tensor de inércia do corpo rígido em relação a um referencial arbitrário
com origem no ponto , e representa a velocidade angular do corpo que é coincidente com a
velocidade angular do referencial solidário com o corpo, Gxyz. Neste momento é importante a
escolha do referencial em que se representa o momento angular, e quanto a isso tem-se
essencialmente duas alternativas: ou se escolhe um referencial fixo que permaneça imóvel enquanto
o corpo roda em torno de , ou se usa um referencial que se mova solidariamente com o corpo
rígido, participando em todos os seus movimentos. Atendendo a que vai ser necessário derivar o
momento angular em ordem ao tempo, a segunda alternativa é preferível, visto que então todas as
componentes do tensor de inércia permanecem constantes ao longo do tempo. No entanto, será
necessário ter em conta que o momento angular está escrito num referencial animado de movimento
de rotação, aplicando-se-lhe então as correspondentes regras específicas da derivação em ordem ao
tempo num referencial em rotação (2.9), tem-se então,
(2.20)
sendo a velocidade angular do referencial.
Se além de se escolher um referencial solidário com o corpo, os seus eixos forem
coincidentes com as direcções principais de inércia, a escrita de (2.19) simplifica-se para
(2.21)
assim, o primeiro termo do segundo membro de (2.20) fica
(2.22)
onde , , são as componentes da aceleração angular do corpo rígido no referencial das
direcções principais de inércia em relação ao ponto . Por outro lado, se o referencial se move
11
solidariamente com o corpo a sua velocidade angular coincide com a velocidade angular do corpo ,
sendo assim, o segundo termo do segundo membro de (2.20) escreve-se:
(2.23)
Finalmente, introduzindo estes resultados em (2.17), obtém-se as equações que modelam o
movimento de rotação de um corpo rígido em torno de um ponto fixo. Estas equações são
conhecidas como Equações de Euler [17].
(2.24)
Discretizando as derivadas temporais com o método de Euler explícito de primeira ordem,
nos instantes e ( , para um intervalo de tempo constante, através da relação:
(2.25)
e substituindo-a nas Equações de Euler (2.24),
(2.26)
obtém-se um conjunto de equações que permitem estimar o vector velocidade angular do corpo
rígido no instante ( , a partir do seu valor no instante e do momento resultante que actua no
corpo no instante . Na secção seguinte veremos como é possível, partindo da velocidade angular
do corpo, determinar a sua nova orientação relativamente ao referencial de orientação fixa .
12
2.4. Matriz de Rotação e Ângulos de Euler
Uma matriz de rotação é uma matriz que quando multiplicada por um vector tem como efeito
mudar a sua orientação sem alterar a sua magnitude. Na representação da rotação por ângulos de
Euler considera-se uma matriz de rotação , que pode ser obtida pelo produto de três matrizes de
rotação independentes, correspondentes a três rotações elementares em torno de cada um dos eixos
do sistema de coordenadas Cartesiano. Considerando que essas matrizes associadas a cada um
dos eixos sejam dadas por , e , os ângulos , e , designados por roll ou
spin, pitch e yaw respectivamente, são chamados Ângulos de Euler. Esta é a versão aeronáutica dos
Ângulos de Euler, pois trocando a ordem das rotações obtêm-se sequências de Ângulos de Euler
diferentes, às quais correspondem matrizes de rotação também diferentes associadas a outra ordem
de rotações elementares [18]. Existem ao todo, doze possibilidades diferentes de os definir, no
entanto apenas três destas possibilidades se usam frequentemente nas diferentes áreas de
aplicação. Na presente modelação usou-se a supracitada versão aeronáutica dos Ângulos de Euler.
Figura 3 - Ângulos de Euler
13
Consideremos uma rotação constituída por três rotações elementares, em que inicialmente o
referencial se encontra na posição xyz, vide Figura 3, assim a primeira rotação do referencial é em
torno do eixo dos xx segundo um ângulo de amplitude a qual deixa o referencial na posição x'y'z'.
Após essa rotação o plano móvel x'y' intercepta o plano fixo yz segundo o eixo y'y', que define a
denominada linha dos nós. A segunda rotação faz-se em torno de y'y', segundo um ângulo ,
passando o referencial da posição x'y'z' para a posição x''y''z''. Por ultimo roda-se em torno de z''z''
segundo um ângulo , ficando o referencial na posição x'''y'''z'''. As matrizes correspondentes às três
rotações elementares são então:
(2.27)
(2.28)
(2.29)
Como o produto de matrizes não é comutativo a ordem das rotações intermediárias, no caso
de rotações sucessivas afecta a matriz resultante. Assim, assumindo uma rotação inicial em x, uma
segunda em y e uma rotação final em torno de z, a matriz de rotação resultante é dada por:
(2.30)
(2.31)
O produto desta matriz por qualquer vector que defina uma posição no espaço, roda-o segundo uma
rotação, que corresponde às três rotações elementares dos três Ângulos de Euler pretendidas. Esta
matriz foi usada para fazer a transformação de coordenadas entre o referencial móvel Gxyz e o
referencial de orientação fixa GX'Y'Z', vide Figura 3. Esta matriz foi codificada em linguagem
FORTRAN, numa rotina externa ao programa Star-CD usado na execução desta Tese, como será
descrito mais detalhadamente adiante.
14
2.5. Rotação do Corpo Rígido – Derivada da Matriz Rotação
Da definição do produto externo entre dois vectores, é possível definir o operador matricial
anti-simétrico . Este operador é útil para efectuar o produto externo da velocidade angular por
outro vector através da multiplicação de uma matriz por esse vector.
(2.32)
(2.33)
Portanto, definindo
(2.34)
(2.35)
Partindo do princípio que se conhece a matriz de rotação que permite rodar um vector escrito
no referencial de orientação fixa GX'Y'Z', para o referencial móvel Gxyz e que esta tem nas suas
colunas os versores dos eixos do referencial móvel escritos nas coordenadas do referencial fixo,
pode-se aplicando (2.9), calcular a derivada temporal destes versores, uma vez que estes apenas
rodam, não variando de qualquer outra forma no tempo. Separando em colunas da matriz
(2.36)
(2.37)
e aplicando o operador (2.34) pode-se escrever o produto externo como um produto matricial
15
(2.38)
aplicando este operador a todas as colunas de obtém-se a derivada temporal da matriz rotação.
(2.39)
(2.40)
Uma vez obtido este resultado, considerando e (com ) as entradas das matrizes e
respectivamente, podem-se discretizar as derivadas , utilizando o método de Euler explícito de
primeira ordem da seguinte forma:
(2.41)
(2.42)
As relações, (2.40) e (2.42), são as que vão permitir estimar a matriz de rotação que
determina a orientação do corpo no instante numérico seguinte. A técnica que se tem vindo a expor é
apenas uma das formas conhecidas de modelar computacionalmente o movimento tridimensional de
um corpo rígido, outra via de abordagem ao assunto seria através da utilização de quaterniões, o que
não foi objectivo deste trabalho. O uso de matrizes de rotação e de quaterniões constituem as
técnicas mais usadas em software no que diz respeito ao movimento tridimensional de objectos,
como é o caso particular dos videojogos. Em caso de haver interesse em estes ou outros assuntos
relacionados com a simulação computacional da dinâmica do corpo rígido, estes podem ser
aprofundados por exemplo em [19,20].
16
3. Modelação do Escoamento
3.1. Equações Base
Em virtude deste estudo incidir sobre escoamentos com velocidades pequenas (sempre
inferiores a 1m/s), o número de Mach será também pequeno, nesse caso pode-se desprezar
variações da densidade do fluido e o escoamento considera-se incompressível. Por outro lado,
considera-se que o fluido é Newtoniano, ou seja, que as tensões viscosas são proporcionais às taxas
de deformação e à viscosidade. As equações que modelam os princípios fundamentais da mecânica
dos meios contínuos, são a da conservação de massa e as de balanço de quantidade de movimento.
Estas últimas, com as simplificações referidas acima têm a designação de equações de Navier-
Stokes, que em notação indicial se escrevem:
(3.1)
(3.2)
em que, os vários símbolos que nelas figuram representam:
Campo de velocidades em cada instante
Pressão estática relativa à hidrostática local
Densidade do fluido
Viscosidade cinemática do fluido,
Nesta modelação o número de Reynolds terá valores próximos de 1200, segundo a definição
(5.3), que em geral maximiza o número de Reynolds local. Note-se que quando o corpo ataca
longitudinalmente o escoamento o comprimento característico deveria ser a espessura (vide secção
5.2 e Figura 36). E quando o corpo cai transversalmente é quando vai com velocidade mais baixa.
Assim estamos perante um escoamento transiente com números de Reynolds locais várias vezes
inferiores ao número de Reynolds nominal. Desta forma o escoamento para efeitos da presente
17
análise, é tido como laminar. Esta abordagem está em conformidade com os trabalhos de referência
usados nesta Tese para efeitos de comparação de resultados [8,9,10,12 ].
3.2 Método de Discretização
As equações diferenciais que modelam o escoamento, já vistas anteriormente, ou seja as
equações de conservação de massa e de balanço de quantidade de movimento, são discretizadas
com o método do volume finito (MVF), conforme implementação no código do Star-CD [21].
Historicamente, o método do volume finito foi introduzido no campo da dinâmica dos fluidos
computacional em 1971 na forma bidimensional e em 1973 na forma tridimensional. Enquanto que o
método das diferenças finitas se baseia na substituição dos operadores diferenciais por outro
algébricos obtidos a partir da série de Taylor, o MVF tem por base a integração das equações de
conservação de massa e balanço de quantidade de movimento. Este método é conservativo por
construção, ou seja, baseia-se na conservação das entidades físicas em cada volume de controlo,
portanto também é conservativo na globalidade. O MVF é particularmente útil quando o domínio
computacional é discretizado por malhas não uniformes, especialmente em situações bi e
tridimensionais, sendo portanto mais adequado para geometrias complexas [22,23].
3.2.1 Esquemas de Discretização do Termo Convectivo
Esquemas de primeira ordem: Estes esquemas geram equações que são mais leves
computacionalmente do que os esquemas de ordem superior, dada a sua simplicidade matemática,
tendo contudo um efeito indesejável denominado de difusão numérica. Este efeito pode ser diluído
com a refinação da malha, o que por outra via também provoca um cálculo mais pesado. Na prática
tem-se que impor um compromisso entre a complexidade do esquema de discretização e o
refinamento das malhas utilizadas. No entanto, para a mesma malha o uso de esquemas de ordem
superior permite (em princípio) melhores resultados.
Esquemas de alta ordem preservam melhor os gradientes mais intensos, mas originam
sistemas de equações mais pesados para resolver, podendo em alguns casos gerar instabilidade
numérica (divergência) e/ou originar soluções exibindo oscilações espaciais sem realidade física
(dispersão numérica) [23].
18
O programa Star-CD tem implementado no seu código os seguintes esquemas de discretização do
termo convectivo da equação de balanço de quantidade de movimento [21]:
UD (Upwind differencing); Esquema de primeira ordem que utiliza informação das células a
montante.
LUD (Linear Upwind Differencing); Este é um esquema adaptado do UD, mas de segunda
ordem, formulado para malhas não estruturadas. Não origina tanta difusão numérica como o
UD, podendo no entanto originar alguma dispersão numérica na presença de gradientes
intensos.
CD (Central Differencing); Esquema de segunda ordem, mas utiliza informação de células
colocada a montante e a jusante do escoamento.
MARS (Monotone Advection and Reconstruction Scheme); Este é um esquema multi-
dimensional de segunda ordem que opera em dois passos distintos: Reconstrução e
Advecção. De todos os esquemas este é o que mostra menos sensibilidade à estrutura da
malha e sua qualidade.É um esquema de segunda ordem que utiliza um plano interpolado
pelo método dos mínimos quadrados.
QUICK (Quadratic Upstream Interpolation of Convective Kinematics); Esquema de segunda
ordem que utiliza uma função quadrática.
Blended differencing; Esquema baseado numa ponderação entre UD com MARS, tendo o
utilizador acesso ao factor de ponderação usado.
Na realização da presente modelação foi usado principalmente o esquema UD, a não ser na
fase das simulações preliminares, na qual se realizaram algumas simulações usando o MARS. No
entanto, não obtendo resultados significativamente melhores que justificassem o uso do MARS,
optou-se pelo UD, pois que, para obter um tempo razoavelmente curto de cálculo sem descurar a
qualidade dos resultados, há que se estabelecer um compromisso entre o esquema usado e o
refinamento da malha.
19
3.2.2 Método de Acoplamento Pressão Velocidade
Os algoritmos disponíveis para o acoplamento pressão velocidade no Star-CD são [21]:
SIMPLE – (semi-implicit method for pressure-linked equations). Disponível apenas para
problemas estacionários.
SIMPISO – (semi-implicit method for pressure-implicit split-operator). Disponível apenas para
problemas estacionários.
PISO – (pressure-implicit split-operator). Único disponível pelo programa para regime
transiente, mas pode ser usado em regime estacionário.
Estes algoritmos são procedimentos iterativos para resolver o sistema de equações
composto pela equação da continuidade e a de balanço de quantidade de movimento e permitem
calcular a velocidade e a pressão. Em regime transiente usa-se PISO e em regime estacionário
pode-se optar por qualquer um deles. Todos se baseiam em estimar a pressão e a velocidade
através dos valores da etapa anterior e seguidamente corrigindo-as, sendo que a principal diferença
é que SIMPLE apenas faz uma correcção e PISO faz várias, usualmente não mais de quatro. As
estimativas para a pressão e velocidade do algoritmo PISO já em si são melhores que as do SIMPLE
devido ao facto de usarem mais termos dos equações de balanço de quantidade de movimento. O
algoritmo SIMPISO é semelhante ao SIMPLE, mas usa uma aproximação melhor para as estimativas
da pressão e velocidade. No presente trabalho, a não ser para algumas experiências iniciais em
regime estacionário onde se usou SIMPLE, para a modelação efectuada usou-se o algoritmo PISO
por virtude de se estar em regime transiente.
3.2.3 Esquemas de Discretização do Termo Temporal
Quanto à dicretização temporal, o Star-CD disponibiliza duas opções: implícito e Crank-
Nicholson. O esquemas implícito para o avanço no tempo, é de primeira ordem, denominado
esquema de Euler totalmente implícito. Esta aproximação usa a suposição da variação linear da
variável entre os dois níveis de tempo. Sendo uma formulação totalmente implícita permite maiores
passos temporais que formulações explícitas, o método está ilustrado na Figura 4. Quanto ao
esquema de Crank-Nicholson é um método de segunda ordem semi-implícito que é usado para
resolver equações diferencias às derivadas parciais [21]. Na presente Tese utilizou-se o esquema de
Euler totalmente implícito.
20
Figura 4 - Molécula Computacional para o Esquema Implícito (Fonte: [21])
Nota: O termo difusivo da equação de balanço de quantidade de movimento é discretizado
pelo método das diferenças centrais e o programa usado apenas tem este esquema implementado
para a discretização deste termo.
21
3.3. Geração de Malha
3.3.1. Condições de Fronteira Utilizadas
O programa usado nesta Tese, Star-CD, tem disponíveis a maioria das condições de
fronteira (boundary conditions) necessárias nas aplicações usuais. As condições de fronteira são
condições que se impõem nas regiões fronteira, estas são constituídas por um conjunto de
elementos bidimensionais especiais, que consistem em algumas faces dos elementos de volume
tridimensionais que formam a malha geral. São formadas por uma superfície (shell) que envolve a
parte do volume pretendido. As condições de fronteira suportados pelo Star-CD que foram usadas
são:
wall Parede, condição que como o nome indica, permite simular uma parede. No
Star-CD, é possivel indicar se a parede é permeável e se é de escorregamento. Neste
trabalho apenas se usaram paredes impermeáveis e de não escorregamento.
inlet Entrada, indicada para as fronteiras onde o fluido entra no domínio
computacional, usada principalmente quando se conhece o perfil de velocidade na entrada.
outlet Saída, região onde os gradientes de velocidade ao longo da direcção do
escoamento são considerados nulos e condicionadas pelas condições de continuidade em
termos globais.
pressure Pressão constante, condição de fronteira em que se considera que a pressão
piezometrica é constante, condição de fronteira amplamente utilizada em regiões fronteira de
saida. Nesta Tese foi utilizada em algumas fronteiras onde se previa que o fluido
abandonasse o domínio simulado.
symmetry Plano de simetria, é equivalente a colocar uma parede de escorregamento, o
mesmo que uma parede de reflexão. O código do Star-CD, é um código intrinsecamente
tridimensional, este tipo de condição fronteira é útil para simular situações bidimensionais.
Este tipo de condição foi usado em todas as modelações bidimensionais apresentadas nesta
Tese.
22
attach Ligação deslizante, consiste numa ligação entre partes separadas da malha,
que permite a livre circulação de fluido entre elas, mas permitindo também que estas regiões
de fronteira deslizem entre si. Usam-se sempre aos pares, uma de cada lado da superfície
deslizante. Nesta Tese este tipo de fronteira foi amplamente utilizada, uma vez que a zona
mais interior do domínio computacional é móvel, movendo-se solidariamente com o objecto
modelado, situado no seu interior.
3.3.2 Mecanismos de Geração de Malha
Essencialmente, existem três mecanismos para crias malhas no programa Star-CD:
Podem-se criar malhas através da linha de comandos do programa introduzindo instruções
nativas do Star-CD e instruções script TCL no módulo ProStar. Tal processo consiste em
criar vértices e linhas que unem os vértices, depois superfícies que ficam limitadas por essas
linhas então discretizam-se as superficies e por fim discretizam-se os volumes por entre
essas superfícies. Uma vez criada a malha pode-se criar também através de comandos
específicos, regiões de fronteira e definir condições de fronteira. Finalmente o ProStar possui
ainda comandos que permitem verificar se a malha está bem definida. É um processo
moroso e pouco prático para malhas complexas.
Outro método consiste em elaborar um script, que é essencialmente como o processo
anterior, usando os mesmos tipos de instruções, mas registando todas elas sequencialmente
para as fazer correr em blocos ou todas de uma vez. Este processo foi o usado para fazer a
malha estruturada que esta modelação usa nas simulações bidimensionais.
Por último e indicado para geometrias complexas, geralmente tridimensionais, usa-se o
processo ilustrado no diagrama da Figura 5. Este foi o método usado para elaborar a malha
usada na simulação 3D. Consiste em começar por modelar a geometria num programa
denominado por StarDesign, onde se modelam sólidos, no presente trabalho, um disco,
esferas e um prisma quadrangular. Depois passa-se para outro programa, o ProSurf, onde se
discretizam as superfícies, seguidamente usa-se o Pro-AM no qual se cria a malha nos
volumes e por fim o ProStar. As regiões de fronteira tanto se podem criar no Pro-AM como
23
no ProStar, no entanto as condições de fronteira e as outras características do problema
definem-se no ProStar. Esta foi a técnica usada para fazer a malha da modelação
tridimensional. No presente caso exigiu um cuidado acrescido devido à existência de uma
superfície de deslizamento, attach, que consiste em duas superfícies esféricas sobrepostas
que deslizam uma em relação a outra. São estas superfícies que unem as duas malhas, a
que tem movimento de rotação e translação (interior) e a que apenas tem movimento de
translação (exterior), vide Figuras 15 e 16.
STAR
PROSTAR
STAR DESIGN
PROSURF
PRO-AM
RotinasExternas
Motor de Cálculo
{
{
{
{
- Modelo Sólido- Pode receber modelos externos
- Tratamento de superfície- Discretização de superfície
- Pode receber superfícies do exterior
- Gerador automático de malha em volume- Pode definir fronteiras- Pode receber malhas de superfície do exterior
- Pré e pós processamento- Capacidade de gerar malha
Figura 5 - Método de Geração de Malha Automática
24
3.3.3. Malha 2D
No caso de se pretender estudar um escoamento exterior, a primeira condição para definir as
dimensões do domínio computacional é que a sua fronteira esteja suficientemente longe para não
interferir com o fenómeno a estudar, neste caso decidiu-se que a dimensão menor do domínio
computacional seria 45 vezes superior à maior dimensão do corpo a simular.
Figura 6 - A Malha; B Cotagem da Malha
Neste estudo, o objecto que se pretende modelar é uma placa plana bidimensional, com uma
espessura de 0.00081m, e com um comprimento 8 vezes superior. As dimensões desta placa bem
como da malha que a envolve podem ser consultadas no diagrama da Figura 6. Esta malha é
constituída essencialmente por três zonas:
(i) Zona central, a vermelho na Figura 6, é a zona da malha com as células mais refinadas.
Nesta parte da malha as células estão dispostas segundo uma evolução radial de acordo
com uma progressão geométrica de razão 1.02, como se mostra nas Figuras 7, 8 e 9. Uma
25
vez que se trata de uma malha estruturada por blocos foi concebida através da elaboração
de um script, tal como se descreveu na secção anterior 3.3.2, começando por se definir
vértices, unindo-os seguidamente e por fim preenchendo a área criada com células como se
pode observar na Figura 8. Sendo assim, à volta da parede da placa existe uma malha
extremamente fina que aumenta de dimensão também segundo uma progressão geométrica
de dentro para fora, (Figura 9).
Figura 7 - Pormenores da Parte Central da Malha 2D
Figura 8 - Método de Elaboração da Malha 2D
26
Figura 9 - Pormenor das Células Junto à Placa
(ii) Coroa circular a azul nas Figuras 6, 7 e 10. Por fora da zona da malha descrita anteriormente
existe uma parte da malha com forma de coroa circular constituída por células de dimensão
maior que as da zona anterior, estas também aumentam de dimensão ao longo da direcção
radial segundo uma progressão geométrica, mas agora de razão 1.03. O conjunto das
malhas representadas a azul e a vermelho, na Figura 7, constitui a parte móvel da
malha, enquanto que a parte representada a amarelo constitui a parte fixa (Figura 10), isto
no que diz respeito à rotação porque em relação ao movimento de translação toda a
malha se pode mover solidariamente.
(iii) Zona a amarelo das Figuras 6, 7 e 10. Por último e envolvendo as zonas descritas acima,
encontra-se a parte exterior da malha, esta com células bastante maiores que as das zonas
descritas anteriormente. O limite interior desta zona de malha tem forma circular enquanto
que o limite exterior tem forma rectangular. A dimensão das células desta parte da malha
também aumenta segundo uma progressão geométrica com razão 1.04, de dentro para fora
na direcção radial, no que diz respeito ao quadrado central (Figura 10), e na direcção
vertical nos rectângulos superior e inferior evidenciados na Figura 6-A.
27
Figura 10 - Pormenor da Parte Fixa da Malha
Como foi referido anteriormente a criação desta malha foi realizada através de um script, no
qual se começou por definir alguns parâmetros, como por exemplo uma dimensão característica da
placa e o número de células pretendido ao longo da sua espessura, em relação aos quais tudo é
calculado com dimensões relativas, tal procedimento tem por objectivo poder alterar algumas
características da malha com bastante economia de tempo. A malha tal como foi utilizada nas
simulações da presente Tese ficou com 48640 células, das quais 36352 pertencem à parte com
capacidade de rotação, ou seja as zonas (i) e (ii) descritas acima.
Na interface entre a zona azul e amarela existe uma superfície que desliza e permite que a
malha móvel se mova no interior da malha fixa, esta funcionalidade pode ser observada nas Figuras
11 e 12, este deslizamento é realizado através de duas superfícies cilíndricas sobrepostas que se
podem mover uma em relação à outra. Tal é realizado no Star-CD recorrendo ao conceito de evento
implementado neste e ao uso de rotinas externas [24], desenvolvidas no decurso deste trabalho, que
calculam em cada iteração a nova posição dos vértices da parte móvel da malha (Secção 4), é
nestas rotinas que está implementada a mecânica do corpo rígido.
28
Figura 11 - Pormenor da União da Malha Móvel e Malha Fixa
Figura 12 – Várias Posições da Parte da Malha com Capacidade de Rotação
Nesta malha as condições de fronteira utilizadas foram as seguintes:
Superfície de escorregamento (Attach) entre a malha móvel e malha fixa (Figura 13-A), esta
condição de fronteira foi implementada na superfície de fronteira entre as zonas
representadas a azul e a amarelo na Figura 6.
Planos de simetria (symmetry) à frente e atrás da malha para simular efeito bidimensional
(Figura 13-B), em relação à direcção perpendicular ao plano da Figura 6.
29
Condição de entrada (inlet) na superfície do lado direito, na superfície do lado esquerdo e na
superfície inferior em relação ao diagrama da Figura 6. O perfil de velocidades utilizado em
todas as entradas foi constante e igual a zero, seria outro valor no caso de se pretender
simular vento lateral por exemplo. No entanto esta condição de fronteira permite que o fluido
a atravesse por forma a que o deslocamento da malha não altere o escoamento que se
pretende.
Condição de pressão estática constante (pressure) na superfície superior (Figura 13-B).
Condição de parede (wall) na superfície de fronteira da placa (Figura 13-C).
Nota: O significado de cada condição fronteira foi já descrito na Secção 3.3.1.
Figura 13 - Condições Fronteira da Malha 2D: A - Attach; B - Entrada, Saída e Semi-Plano; C - Parede
30
3.3.4. Malha 3D
Figura 14 - Domínio da Malha Tridimensional
Para gerar esta malha tridimensional, ao contrário da malha bidimensional anterior, utilizou-
se software de apoio à geração automática de malha. Primeiro modelaram-se os sólidos com auxílio
do software StarDesign. Seguidamente procedeu-se ao tratamento e discretização das superfícies
dos sólidos no ProSurf, para então se gerar a malha em volume no Pro-AM. Finalmente, no ProStar,
definiram-se condições do problema tais como: regiões de fronteira, condições de fronteira e
propriedades do fluido, tal método está exposto no diagrama da Figura 5.
Figura 15 - Cortes da Malha 3D
31
Esta malha contém no seu interior um disco com 1 centímetro de diâmetro e 1 milímetro de
espessura, o disco está representado a negro nas Figuras 15 e 16, e a branco na Figura 17. A
malha é constituída por camadas de malha cada vez menos refinadas, à medida que nos afastamos
do centro (disco), e contém uma superfície esférica de escorregamento (attach) entre as zonas
assinaladas a azul e a amarelo nas Figuras 15 e 16. O domínio externo tem a forma de um prisma
quadrangular (Figura 14), cujas dimensões são 100 vezes a espessura do disco para a aresta da
base do prisma, e 200 vezes a espessura do disco para a altura, pois é nessa direcção que se prevê
que fique a parte principal da esteira do escoamento causada pelo movimento do disco.
Figura 16 - Parte Interior da Malha
Esta malha é constituída por uma parte móvel, representada a azul e vermelho nas Figura 15
e 16, que tem forma esférica para que possa rodar em torno de qualquer eixo arbitrário. E por uma
parte fixa exterior representada a amarelo e verde nas Figuras 15 e 16. O deslizamento entre elas é
realizado através de duas superfícies esféricas sobrepostas com raio igual a 6 vezes o raio do disco,
que se podem mover uma em relação à outra. Tal é realizado no Star-CD recorrendo ao conceito de
evento implementado neste e ao uso de rotinas externas [24], desenvolvidas no decurso deste
trabalho que calculam em cada iteração a nova posição dos vértices da parte móvel da malha
(Secção 4), é nestas rotinas que está implementada a mecânica do corpo rígido.
32
Figura 17 - Pormenor da Malha em Volta do Disco
Na Figura 17 evidencia-se a dimensão relativa das células usadas junto às paredes do disco,
estas foram construídas de forma a aumentarem de tamanho na direcção exterior segundo uma
progressão geométrica de razão 1.05. A escolha da dimensão das células teve em conta que não se
pode usar células tão pequenas quanto seria desejado para não comprometer o tempo de cálculo e
as limitações dos computadores utilizados. Mesmo assim, para a simulação realizada com esta
malha, para se obter um segundo de simulação usaram-se 6 processadores em paralelo durante 4
dias. Esta malha foi usada com um total de 969320 células das quais 228064 estão na parte com
movimento de rotação.
Nesta malha as condições de fronteira utilizadas foram as seguintes:
Superfície de escorregamento (Attach) entre a malha móvel e malha fixa, ou seja entre as
conchas esféricas concêntricas representadas em corte a azul e a amarelo na Figura 16.
Condição de entrada (inlet) em todas as superfícies laterais do prisma e também na base
inferior (Figura 14). O perfil de velocidades utilizado em todas as entradas foi constante e
igual a zero, seria outro valor no caso de se pretender simular um escoamento de
aproximação diferente. Esta condição de fronteira permite que o fluido a atravesse por forma
a que o deslocamento da malha não altere o escoamento que se pretende.
Condição de pressão estática constante (pressure) na superfície de fronteira da base
superior do prisma.
Condição de parede (wall) na superfície de fronteira constituída pelas paredes do disco
(Figura 17).
Nota: O significado de cada condição fronteira foi já descrito na Secção 3.3.1.
33
4. Algoritmo de Cálculo
Figura 18 - Algoritmo de Cálculo
O Algoritmo de cálculo é na sua essência um ciclo de fases sucessivas, em que cada volta
constitui uma iteração da simulação em execução (Figura 18). O controlo geral do processo é
assegurado pelo programa usado nesta modelação Star-CD, além disso este é responsável pela
parte mais complexa do calculo, que é o cálculo do escoamento descrito sumariamente na secção 3.
Após cada iteração do cálculo do escoamento a execução passa para a rotina posdat na qual se
calculam as forças e os momentos exercidos no corpo, a força é calculada através do somatório das
Cálculo do Escoamento
Cálculo das Forças
Cálculo dos Momentos
Cálculo da Velocidade
de Translação
Cálculo da Posição do Centro de
Massa
Cálculo da Velocidade Angular
Equações de Euler
Cálculo da Matriz de Rotação
Ortonormalização da Matriz de
Rotação
Escrita de Resultados em
Ficheiro
Alteração das Coordenadas
dos Vértices da Malha
34
forças em cada célula bidimensional da superfície do corpo, estes valores são disponibilizados pelo
Star-CD, o momento é calculado de acordo com (4.1).
(4.1)
(4.2)
Posteriormente, já com estes resultados disponíveis numa outra rotina externa newxyz,
(Figura 19), fazem-se todos os cálculos do movimento quanto a translação, essencialmente
recorrendo às relações (2.15) e (2.16), ficando a conhecer-se a velocidade do corpo e a nova
posição do centro de massa. Em seguida e também nesta rotina fazem-se os cálculos relativos à
rotação do corpo através das Equações de Euler (2.26), calcula-se a velocidade angular, e também a
nova matriz de rotação de acordo com (2.42). Ainda nesta rotina (newxyz) procede-se à
ortonormalização da matriz de rotação e ao armazenamento dos resultados em ficheiros de dados.
Por fim esta rotina altera as coordenadas dos vértices da malha, começando por alterar as
coordenadas dos vértices da parte que roda e posteriormente de toda a malha no que diz respeito à
translação, voltando o programa ao cálculo do escoamento, portanto a uma nova iteração.
35
Colocar Corpo naPosição Inicial
Dinâmica da Translação
Dinâmica da Rotação(Equações de Euler)
Alteração das Posiçõesdos Vértices da Malha
Star-CD
Cálculo do Escoamento
Controlo Geraldo Processo
posdat.f
newxyz.f
bcdefw.f bcdefi.f
Cálculo do EscoamentoModelção Tridimensional
da Dinâmica de Corpo Rígido
Cálculo de Forças e Momentos
Figura 19 - Rotinas Externas
As rotinas bcdefw e bcdefi existem com a função de indicar as velocidades do escoamento
junto às regiões de fronteira definidas como paredes e definidas como entradas, respectivamente.
Na primeira iteração do processo descrito acima, a rotina newxyz coloca o corpo objecto da
simulação na posição inicial pretendida, recorrendo à matriz de rotação com as entradas descritas
em função dos Ângulos de Euler, como descrito na relação (2.31).
As rotinas descritas acima, representadas na Figura 19, foram codificadas em linguagem
FORTRAN propositadamente para o fim pretendido nesta Tese.
36
5. Resultados
5.1. Simulações Preliminares
Na fase inicial deste trabalho realizaram-se várias simulações de teste para afinar o processo
geral. Apresentam-se aqui apenas algumas notas sobre esses testes preliminares.
Feito o código da rotina em FORTRAN que codifica a dinâmica do corpo rígido, este foi
testado, quer quanto à translação quer quanto à rotação. Estes testes foram realizados ainda sem a
modelação do escoamento, ou seja, estava-se a modelar a dinâmica de corpos sem interacção com
nenhum fluido. Nestes testes foram simulados lançamentos de projécteis variando as condições
iniciais, ângulos e velocidade inicial, e também lançamento de projécteis de acordo com o modelo
designado por point mass [25]. Para verificar a rotina quanto à rotação foi simulado um pêndulo
gravítico a oscilar em vários planos. Foram ainda realizadas outras imposições nas forças e/ou
momentos e também nas velocidades iniciais, linear e angular. Contudo, devido à falta de resultados
teóricos para a maior parte das situações, apenas se pode dizer que se obteve para todas as
simulações resultados fisicamente credíveis. Verificando-se também que ao baixar o passo temporal
o erro relativo dos resultados baixava da mesma forma, o que era de esperar uma vez que o
método é de primeira ordem.
Quando se codificou a dinâmica do corpo rígido, usou-se um passo genérico para o número
de iterações ao fim das quais se ortonormaliza a matriz de rotação. Dos vários testes realizados
concluiu-se que é preferível ortonormalizar a matriz de rotação em todas as iterações.
Noutra fase do trabalho testou-se o modelo já com a dinâmica do corpo rígido acoplada ao
modelo de cálculo de escoamento, sendo que, numa fase inicial se fizeram simulações em que se
impôs o movimento do objecto, tanto de translação como de rotação. Foi nesta fase que se tomaram
decisões necessárias quanto às opções que o Star-CD disponibiliza, de forma a que o cálculo não
ficasse instável e divergisse, tais como escolha da malha, passo temporal, esquema de discretização
do termo convectivo, condições de fronteira e parâmetro de tolerância do erro da ligação das
interfaces móveis.
37
5.2. Resultados das Modelações Bidimensionais (2D)
As modelações que se pretendem fazer nesta secção são bidimensionais e o objecto que se
pretende modelar é uma placa plana, de secção rectangular, cujo comprimento é oito vezes a sua
espessura. Estas modelações foram todas realizadas com a malha descrita na secção 3.3.3.
Nas simulações bidimensionais apresentadas nesta Tese, sempre que não for dito nada em
contrário, foram usados os parâmetros físicos seguintes nas diversas simulações, que são iguais aos
da referência [12], que serviu de base de comparação para este trabalho. O critério de convergência
entre cada passo temporal usados nestas simulações foi 0.1% do resíduo, que é o valor que o Star-
CD usa por defeito.
Grandeza Símbolo Valor Unidades
Comprimento da placa 0.00648
Espessura da placa 0.00081
Passo temporal 1/4000
Massa específica do fluido 1000
Massa específica do corpo 2700
Viscosidade dinâmica 0.0008887
Aceleração da gravidade 9.8
Magnitude da velocidade inicial 0.12615
Ângulo inicial em relação à
horizontal 45 [grau]
Tabela 1 - Valores dos Parâmetros Físicos Usados por Defeito em Todas as Simulações 2D
Nota: A velocidade inicial é aplicada na placa na sua direcção longitudinal e no sentido descendente.
38
5.2.1 Parâmetros Adimensionais Utilizados na Apresentação dos Resultados
Em sintonia com o Teorema dos Pi’s de Buckingham, pode-se supôr que numa primeira
aproximação, o fenómeno físico em estudo pode ser descrito pelos três números adimensionais
seguintes, isto supondo que o corpo é abandonado com velocidade inicial nula e sempre com o
mesmo ângulo relativamente à horizontal.
(5.1)
(5.2)
(5.3)
(5.4)
onde os símbolos representam:
Razão de densidades
Razão entre a comprimento e a espessura da placa
Número de Reynolds definido em conformidade com a literatura sobre o tema em
questão [8,9,10,12]
Massa específica do corpo sólido
Massa específica do fluido
Comprimento da placa
Espessura da placa
Aceleração da gravidade
Velocidade terminal de queda livre, estimada através do balanço da força gravítica
contra a resistência do fluido numa placa com um coeficiente de resistência unitário [12]
Neste estudo o número adimensional vai manter-se constante, pois que, para o fazer variar
seriam necessárias várias malhas. No entanto decidiu-se manter uma análise semelhante à das
referências sobre o assunto e também para uma mais fácil comparação dos resultados com outros
trabalhos nesta área.
39
Os seguintes números adimensionais, que se podem obter dos resultados das diversas
simulações servirão essencialmente para apresentar os resultados de forma adimensional.
(5.5)
(5.6)
(5.7)
(5.8)
(5.9)
(5.10)
(5.11)
(5.12)
(5.13)
onde os símbolos representam:
Componente da velocidade adimensional na direcção
Componente da velocidade na direcção
Magnitude da velocidade adimensional
Magnitude da velocidade
Componente da velocidade média adimensional na direcção
Velocidade angular adimensional
Velocidade angular
Tempo
Tempo adimensional
Força na direcção
Força adimensional na direcção
Diferença entre a massa do corpo e a massa do fluido com o volume do corpo
Pressão estática adimensionalizada pela pressão dinâmica calculada tendo por base
a velocidade terminal
Pressão estática
40
Posição adimensionalizada na direcção do eixo dos xx
Posição na direcção do eixo dos xx
Posição adimensionalizada na direcção do eixo dos yy
Posição na direcção do eixo dos yy
Comprimento da placa
Com o objectivo de verificar a consistência da análise dimensional foram realizados três
simulações, onde se variaram alguns parâmetros físicos, mas mantendo o valor dos parâmetros
adimensionais, vide Tabela 2. Tendo a simulação A como base de comparação, na simulação B
variou-se apenas a aceleração da gravidade e a viscosidade dinâmica . E na simulação C
alterou-se a aceleração da gravidade e as massas especificas do corpo e do fluido e .
A B C
Valor das grandezas
físicas impostos nas
simulações
Valores
adimensionais
Resultados
adimensionais
Tabela 2 - Verificação da Consistência da Análise Dimensional
41
5.2.2. Validação
Com o intuito de validar o processo foi realizada uma modelação tão próxima quanto possível
da realizada por Changqiu Jin e Kun Xu, em 2008 [12], e dos trabalhos experimentais de Z. Jane
Wang, [8].
Parâmetros físicos da simulação:
Ângulo inicial,
Parâmetros adimensionais:
Figura 20 – Comparação de Trajectórias
42
Como se pode observar na Figura 20, onde estão sobrepostas as trajectórias da presente
modelação e dos trabalhos das referências supracitadas, estas estão relativamente próximas umas
das outras, e estão apresentadas em coordenadas adimensionalizadas para uma mais fácil
comparação com outros trabalhos. Note-se também que nos gráficos das Figuras 21 e 22, o tempo
não começa em zero uma vez que se apresentam comparações de resultados e as referências não
tinham disponíveis resultados para os instantes iniciais, pois o que se pretende mostrar é a
comparação das forças instantâneas na fase do movimento em que este já está estabilizado, fase
terminal. A placa faz o movimento conhecido como tumbling com período bastante semelhante entre
os vários resultados. Também se observa visualmente no gráfico que o declive médio da descida é
aproximadamente igual. Repare-se que para os últimos períodos simulados as discrepâncias são
maiores, o que é natural que aconteça devido à acumulação temporal dos erros do cálculo.
Seguidamente, nas Figuras 21 e 22, apresentam-se as componentes das forças instantâneas
adimensionalizadas que actuam na placa, as do presente cálculo e as da referência [12]. Como se
pode observar nos gráficos a semelhança é muito significativa.
Figura 21 – Força na Direcção y
43
Figura 22 – Força na Direcção x
Nos gráficos das Figuras 23 e 24, apresentam-se os resultados adimensionalizados das
componentes horizontal e vertical da velocidade instantânea e velocidade angular instantânea. Para
estes resultados os trabalhos de referência desta Tese não têm valores disponíveis para
comparação.
44
Figura 23 - Velocidades
Figura 24 – Velocidade Angular
Vx*Vy*
45
5.2.3. Relação Entre Velocidades Terminais e Número de Reynolds.
Nesta secção apresentam-se os resultados de uma série de simulações que se realizaram
variando os parâmetros físicos aceleração da gravidade e viscosidade do fluido, com o objectivo de
variar o número de Reynolds, de forma a manter fixos , B=2.7 (Figura 25) e B=5 (Figura 26).
Com os resultados obtidos calcularam-se: a magnitude da velocidade média; a componente vertical
da velocidade média e a velocidade angular. Os gráficos das Figuras 25 e 26 apresentam as
relações das velocidades descritas com o número de Reynolds definido como em (5.3), para dois
valores distintos da razão de massas especificas B=2.7 e B=5.
Figura 25 - Resultados para e
0
0,2
0,4
0,6
0,8
1
1,2
1,4
0 500 1000 1500 2000 2500
Re*
Vy*
Vm*
Wz*
46
Figura 26 – Resultados para e
Para algumas simulações existem apenas valores para o Vy*, uma vez que uma velocidade
média vertical é sempre possível de calcular, mesmo quando a trajectória é aparentemente caótica.
Os pontos indicados no gráfico da Figura 25, parecem indicar alguma tendência perceptível,
contudo não muito evidente.
No gráfico da Figura 27 fez-se o raciocínio contrário, ou seja, fixou-se o Re e variou-se a
razão de densidades B. Não se conseguiu obter muitos resultados uma vez que quando se aumenta
B significativamente, o corpo precisa de muito tempo para acelerar e atingir uma velocidade próxima
da sua velocidade terminal e até isso acontecer a sua trajectória apresenta aparência caótica.
0
0,2
0,4
0,6
0,8
1
1,2
1,4
0 500 1000 1500 2000 2500 3000 3500
Re*
Vy*
Vm*
Wz*
47
Figura 27 - Resultados para α=8 e =1198
Quando se pretende simulações com relações de densidades B menores, o algoritmo fica
instável, e exige maior poder computacional. E quando se obtêm resultados verifica-se que a
trajectória é frequentemente caótica o que já tinha sido observado em outros estudos sobre o tema
presente [26].
0
0,02
0,04
0,06
0,08
0,1
0,12
2 2,5 3 3,5 4 4,5 5
Vy*
B
48
5.2.4. Simulação com Movimento de Fluttering
Nesta secção mostram-se alguns resultados de uma simulação cuja trajectória têm
características do movimento do tipo fluttering, estes resultados foram obtidos para B=1.89 e
Re*=102. Note-se que para valores de B abaixo de 2 o movimento estabiliza com trajectórias do tipo
fluttering, mas por vezes é necessário muito mais tempo de simulação para obter regime estável, e
para valores de B inferiores a 1.7 o algoritmo utilizado fica instavel.
Figura 28 - Trajectória do Movimento Fluttering
Observe-se que a trajectória mostrada na Figura 28 é semelhante á apresentada na Figura
1(a), contudo com oscilações maiores, movimento característico do fluttering. Nas Figuras 29 e 30
apresentam-se os resultados para as componentes das forças adimensionalizadas, e da velocidade
angular adimensionalizada, respectivamente, ambas em função do tempo adimensional.
49
Figura 29 - Componentes das Forças Adimensionalizadas
Figura 30 - Velocidade Angular Adimensional
5.2.5. Estudo das Trajectórias em Função da Razão de Densidades
Na Figura 31 apresentam-se os resultados de um estudo feito sobre a evolução das
trajectórias da placa em função da razão de densidades B, mantendo-se tudo o resto constante e
com os valores apresentados na Tabela 1. O gráfico da Figura 32 mostra o pormenor do início das
trajectórias, sendo os resultados os mesmos dos da Figura 31.
50
Na Figura 31, a trajectória para B=2.7, é precisamente a da simulação usada para validar
este trabalho, em 5.2.2., simulação tida como padrão de comparação nesta Tese. Observe-se que à
medida que se sobe o valor de B a partir de B=2.7, até cerca de B=5, as trajectórias entram num
regime estável mais rapidamente. E para valores de B maiores a placa inicialmente desce muito
rapidamente, provavelmente até ganhar velocidade suficiente e consequente sustentação, entrando
depois num regime mais estável com rotação rápida. Note-se também que todas as simulações
destas figuras foram feitas para um segundo de tempo total, e as trajectórias para B maiores são
mais longas o que permite observar que têm velocidades mais elevadas.
Figura 31 - Trajectórias em Função de B
51
Figura 32 - Pormenor Inicial das Trajectórias em Função de B
5.2.6. Estudo das Trajectórias em Função do Ângulo de Lançamento Inicial
Na Figura 33 apresentam-se os resultados de um estudo sobre a evolução das trajectórias
da placa em função do ângulo de lançamento inicial , mantendo-se tudo o resto constante e com
os valores descritos na Tabela 1. O gráfico da Figura 34 mostra o pormenor do início das trajectórias,
sendo as simulações as mesmas que estão apresentadas na Figura 33.
Tal como no estudo anterior, salienta-se que a trajectória para 0.25pi é a da simulação que
foi usada na validação desta Tese, mantendo-se também aqui como referência. É importante não
esquecer que os restantes parâmetros destas simulações são os da Tabela 1, em particular que a
velocidade inicial não é zero e é aplicada na direcção longitudinal da placa, sentido descendente,
portanto para a esquerda.
52
Note-se que para as trajectórias com de 0.05pi a 0.20pi, a placa dá uma volta sobre si
própria relativamente cedo e estabilizam para a direita. E que nas trajectórias para de 0.25pi a
0.45pi a placa dá a primeira volta sobre si própria mais tarde e estabilizam para o lado esquerdo.
Figura 33 - Trajectórias em Função do Ângulo Inicial
53
Figura 34 - Pormenor das Trajectórias em Função do Ângulo Inicial
Por fim na Figura 35, apresentam-se algumas trajectórias da simulação com vários ângulos
de lançamento iniciais, mas desta vez com velocidade inicial nula. Observe-se em particular as
trajectórias correspondentes aos 0.00pi e 0.50pi, enquanto que no primeiro caso a placa avança para
o escoamento transversalmente tendo o máximo de resistência (neste caso a resistência é que dá
"sustentação"), no segundo caso a placa avança longitudinalmente para o escoamento tendo o
mínimo de resistência. Repare-se na diferença de comprimentos das trajectórias para o mesmo
intervalo de tempo simulado.
54
Figura 35 - Trajectórias com Velocidade Inicial Nula
55
5.2.7. Campos de Velocidade, de Pressão Estática e Vorticidade
Nesta secção os resultados apresentados são a sequência de campos de velocidades,
pressão estática e vorticidade adimensionalizados. A simulação em questão corresponde à
trajectória da Figura 33 para um ângulo inicial de 0.20pi, trajectória essa com movimento do tipo
tumbling. Os oito instantes de cada uma das Figuras 36, 37 e 38, correspondem ao mesmo instante
pela posição que ocupam em cada sequência. Da primeira imagem à oitava percorrem
aproximadamente um período completo da trajectória em regime estabilizado.
Na Figura 36, apresentam-se os campos de velocidades, note-se que o primeiro e o último
correspondem aproximadamente a um período de diferença temporal.
Na Figura 37, estão apresentados os campos de pressão estática correspondentes à Figura
36, repare-se por exemplo na Figura 37-e, onde se observa alta pressão no "bordo de ataque" da
placa, e baixa pressão no "extradorso".
Por último, na Figura 38, estão campos de vorticidade, observe-se também a semelhança
entre o primeiro e o último e a característica libertação de vórtices associada a uma perda de
sustentação também observada noutros estudos.
56
Figura 36 - Campos de Velocidades Adimensionais. Da esquerda para a direita e de cima para baixo, os instantes (0.44; 0.47; 0.50; 0.53; 0.56; 0.59; 0.62; 0.65) s
57
Figura 37- Campos de Pressões Adimensionais. Da esquerda para a direita e de cima para baixo, os instantes (0.44; 0.47; 0.50; 0.53; 0.56; 0.59; 0.62; 0.65) s
58
Figura 38 - Campos de Vorticidade. Da esquerda para a direita e de cima para baixo, os instantes, (0.44; 0.47; 0.50; 0.53; 0.56; 0.59; 0.62; 0.65) s
59
5.2.8. Outros Resultados Ilustrativos de Simulações 2D
Nesta secção apresentam-se algum resultados que não se enquadravam nas secções
anteriores, na Figura 39 pode-se observar uma ampliação de um campo de velocidades e na Figura
40, o pormenor dos vectores velocidade junto da placa.
Figura 39 - Pormenor de um Campo de Velocidades
Figura 40 - Vectores Velocidade Junto à Placa
60
Figura 41 - Alguns Exemplos de Trajectórias. Da esquerda para a direita e de cima para baixo a, b, c e d
Na Figura 41 apresentam-se alguns tipos de movimentos menos comuns, obtidos para
valores de alguns parâmetros tidos como "menos normais". Figura 41-a, simulação com
. Figura 41-b, simulação obtida com . Figura 41-c, simulação obtida com
.Figura 41-d resultado obtido com e com . Os restantes
valores são os mostrados na Tabela 1.
61
5.3. Resultados da Simulação Tridimensional (3D)
Nesta secção apresentam-se os resultados de uma simulação tridimensional. As
características da malha utilizada estão descritas na secção 3.3.4. Esta simulação exigiu muito mais
poder computacional do que as simulações bidimensionais da secção anterior, por isso apenas se
tem disponível um segundo de simulação. Para realiza-la foi necessário a técnica do processamento
em paralelo, neste caso com seis processadores durante cerca de quatro dias. Não se faz
comparação com outros resultados porque não são conhecidas até então trabalhos experimentais ou
numéricos com os quais se pudesse comparar os resultados.
Parâmetros físicos da simulação:
Ângulo inicial,
Nota: O ângulo inicial de inclinação foi medido em relação ao eixo dos xx, mantendo assim
alguma semelhança com a modelação da placa rectangular nas simulações bidimensionais.
A velocidade inicial aplicou-se no sentido descendente e na direcção correspondente à recta
de intersecção do plano do disco com o plano xoy, (Figura 42).
62
Figura 42 - Referencial Utilizado. Posição Inicial do Disco.
Na Figura 43, apresentam-se duas visualizações das trajectórias tridimensionais do
disco, enquanto que nas Figuras 44 e 45 estão apresentadas as projecções ortogonais
dessas trajectórias.
Figura 43 - Trajectória 3D. Nota: estes eixos não estão normalizados
63
Figura 44 - Trajectórias nos Planos xoy e yoz. Nota: Chama-se a atenção para o facto de a escala não ser a mesma nos dois eixos.
Figura 45 - Trajectórias nos Planos xoz. Nota: Chama-se a atenção para o facto de a escala não ser a mesma nos dois eixos.
64
Seguidamente, nas Figuras 46 a 50, apresentam-se de forma gráfica os seguintes
resultados: as componentes das forças instantâneas que actuam no disco; as componentes dos
momentos instantâneos que actuam no disco; as componentes das velocidades instantâneas e as
componentes da velocidade angular e ainda a variação dos Ângulos de Euler no tempo. Da análise
destes gráficos conclui-se que os resultados se assemelham bastante com a realidade física,
bastando para isso largar um pequeno disco de alumínio num tanque com água, o movimento que se
observa condiz precisamente com os resultados obtidos.
Figura 46 - Componentes da Forças Instantânea que Actua no Disco
65
Figura 47 - Componentes dos Momentos Instantâneos que Actuam no Disco
Figura 48 - Componentes da Velocidade do Disco
66
Figura 49 - Velocidade Angular do Disco
Figura 50 - Ângulos de Euler do Disco
Nas imagens da Figura 51, pode-se observar uma sequência de oito campos de velocidade
em torno do disco em vários instantes igualmente espaçados no tempo, todos eles são cortes do
campo do escoamento do disco pelo plano z=0. Os campos de pressão estática para os mesmos
instantes e segundo o mesmo corte podem ser observados na Figura 52, e os campos de vorticidade
também nos mesmo instantes e corte encontram-se representados na Figura 53. Da análise destas
figuras e das anteriores observa-se que os resultados se assemelham bastante com a realidade
física, ainda que apenas o possamos afirmar de forma qualitativa, bastando para tal abandonar um
pequeno disco de alumínio dentro de um tanque com água.
67
Figura 51 - Campos de Velocidades em m/s. Da esquerda para a direita e de cima para baixo, os instantes (0.71; 0.74; 0.77; 0.80; 0.83; 0.86; 0.89; 0.92) s
68
Figura 52 - Campos de Pressão em Pa. Da esquerda para a direita e de cima para baixo, os instantes (0.71; 0.74; 0.77; 0.80; 0.83; 0.86; 0.89; 0.92) s
69
Figura 53 - Campos de Vorticidade. Da esquerda para a direita e de cima para baixo, os instantes (0.71; 0.74; 0.77; 0.80; 0.83; 0.86; 0.89; 0.92) s
70
Para finalizar apresenta-se ainda em pormenor um campo de velocidades em redor do disco,
corte no plano z=0, aos 0.56 segundos da simulação (Figura 54). E uma representação dos vectores
velocidade do fluido nas imediações do disco ao longo do corte z=0, aos 0.80 segundos da
simulação (Figura 55).
Figura 54 - Pormenor de um Campo de Velocidades
Figura 55 - Vectores Velocidade
71
6. Conclusões
Dos resultados da presente Tese pode-se concluir que é possível realizar o acoplamento entre a
dinâmica do corpo rígido tridimensional e a mecânica dos fluidos computacional, em particular com o
programa comercial Star-CD.
Pode-se modelar o movimento do corpo rígido recorrendo a técnicas de malhas móveis, sendo
que, para isso teve-se que codificar algumas rotinas em FORTRAN especificamente para esse fim.
No entanto esta técnica exige maior poder computacional, que as modelações gerais em que os
corpos sólidos estão fixos ou têm movimentos impostos à partida. Uma vez que as malhas têm que
ser mais refinadas e o passo temporal tem que ser menor. Para resolver esse problema recorreu-
se ao uso de computação em paralelo. Esta técnica foi usada principalmente nas simulações
tridimensionais.
Foi possível identificar em diferentes regimes de escoamento alguns movimentos conhecidos,
como o tumbling e o fluttering, de acordo com a literatura científica sobre o assunto.
O processo usado fornece resultados fisicamente realistas como se pode constatar pelos
resultados obtidos apresentados no Capítulo 5. Foi possível com esta técnica reproduzir os
resultados de alguns estudos recentes sobre queda livre de placas no seio de um fluido, quer
numéricos, quer experimentais [8,12], o que leva a crer que o acoplamento entre a dinâmica do corpo
rígido e a mecânica dos fluidos computacional foi bem sucedido.
72
Referências
[1] Nicolas Sardoy, Jean-Louis Consalvi, Bernard Porterie, A. Carlos Fernandez-Pello, "Modeling
transport and combustion of firebrands from burning trees". Combustion and Flame 150, pp.
151–169 (2007).
[2] Jubaraj Sahu, "Time-Accurate Numerical Prediction of Free Flight", Aerodynamics of
a Finned Projectile, Army Research Laboratory (2005).
[3] S.A. Ansari, Z. Zbikowski, K. Knowles, "Aerodynamic modelling of insect-like flapping flight for
micro air vehicles", Department of Aerospace, Power and Sensors, Cranfield University,
Defence Academy of the United Kingdom, Shrivenham (2006).
[4] David J. Willis, Emily R. Israeli, Per-Olof Persson, Mark Drela, Jaime Peraire, Sharon M.
Swartzk, Kenneth S. Breuer, "A Computational Framework for Fluid Structure. Interaction in
Biologically Inspired Flapping Flight", Applied Aerodynamics Conference (2007).
[5] Nuno André Ramos Maia, "Aerodynamic Characterization of a Micro Air Vehicle, Dissertação
para obtenção do Grau de Mestre em Engenharia Aeroespacial", Intituto Superior Técnico,
(2007).
[6] J. C. Maxwell, "The Scientific Papers of James Clerk Maxwell", Dover, New York, pp. 115-118
(1890).
[7] A. Belmonte, H. Eisenberg, E. Moses, "From flutter to tumble: Inertial drag and Froude
similarity in falling paper", Phys. Rev. Lett., 81 , pp. 345-348 (1998).
[8] A. Anderson, U. Pesavento, Z. Jane Wang, "Unsteady aerodynamics of fluttering and
tumbling plates", J. Fluid Mech., 54, pp. 65–90 (2005).
[9] U. Pesavento, Z.J. Wang, "Falling paper: Navier–Stokes solutions, model of fluid forces, and
center of mass elevation", Phys. Rev. Lett., 93:14 (2004).
[10] A. Anderson, U. Pesavento, Z. Jane Wang, "Analysis of transitions between fluttering,
tumbling and steady descent of falling cards". Journal of Fluid Mechanics, 541 , pp 91-104
(2005).
73
[11] David L. Finn., "Falling Paper and Flying Business Cards", SIAM News, Volume 40, Number 4
(2007).
[12] Changqiu Jin, Kun Xu, "Numerical Study of the Unsteady Aerodynamics of Freely Falling
Plates", Communications in Computational Physics Vol. 3, No. 4, pp. 834-851 (2008).
[13] Daniel Tam, John W.M. Bush., Michael Robitaille, Arshad Kudrolli, "Tumbling dynamics of
flexible wings", Department of Mathematics, Massachusetts Institute of Technology,
Cambridge; and Department of Physics, Clark University, Worcester (2009).
[14] Peter Weiss, "The Puzzle of Flutter and Tumble, Physicists reconsider the fall of leaves",
Science News, vol. 154 No. 18, p. 285 (1998).
[15] Rajat Mittal, Veeraraghavan Seshadri, "Flutter, Tumble and Vortex Induced Autorotation",
Department of Mechanical Engineering, University of Iowa, U.S.A. (2004).
[16] Ferdinand P. Beer, E. Russell Johnston Jr., "Mecânica Vectorial para Engenheiros –
Cinemática e Dinâmica", 5ª Edição. McGraw-Hill (1991).
[17] L. Landau, E. Lifshitz, "Física Teórica: Mecânica", Editora Mir, Moscou (1978).
[18] E. R. Bachmann, R. B. McGhee, X. Yun, M. J. Zyda, "Rigid Body Dynamics, Inertial
Reference Frames and Graphics Coordinate Systems: A Resolution of Conflicting
Conventions and Terminology. For Submission to IEEE Symposium on Computational
Intelligence in Robotics and Automation", Naval Postgraduate School, Monterey, CA 93943
(2006).
[19] David Baraff, "An Introduction to Physically Based Modeling: Rigid Body Simulation I.
Unconstrained Rigid Body Dynamics", Robotics Institute, Carnegie Mellon University (1997).
[20] Ralph Bach, Russell Paielli, "Direct Inversion of Rigid-Body Rotational Dynamics", NASA
Technical Memorandum 102798, (1990).
74
[21] CD-ADAPCO group. Star-CD version 3.26 methodology (2006).
[22] Robert Eymard, Thierry Gallouet, Raphaele Herbin, "Finite Volume Methods" (2006).
[23] Hirsch C., "Numerical Computation of Internal and External Flows, Volume 2: Fundamentals
of Computational Fluid Dynamics", John Wiley & Sons, New York, (2007).
[24] CD-ADAPCO group. Rotating and moving meshes. Star-CD Version 3.26 User Guide, (2005).
[25] Mehmet Akçay, "Development of Universal Flight Trajectory Calculation Method for Unguided
Projectiles", Turkish J. Eng. Env. Sci. 28, pp. 369-376 (2004).
[26] Tomonori Yamada, Shinobu Yoshimura, "Line Search Partitioned Approach for Fluid-structure
Interaction Analysis of Flapping Wing", CMES, vol.24, no.1, pp.51-60 (2008).
[27] White, Frank M., "Mecânica dos Fluidos" McGraw-Hill, (1999).
