Upload
dangthuan
View
229
Download
0
Embed Size (px)
Citation preview
Cálculo da massa adicionada de um corpo com o método
dos painéis
Rui Pedro Silva Caracol
Dissertação para obtenção do Grau de Mestre em
Engenharia Mecânica
Orientador: Prof. José Alberto Caiado Falcão de Campos
Júri
Presidente: Prof. Viriato Sérgio de Almeida Semião
Orientador: Prof. José Alberto Caiado Falcão de Campos
Vogal: Prof. Luís Rego da Cunha Eça
Junho de 2014
i
Agradecimentos
Gostaria de agradecer em primeiro lugar ao professor Falcão de Campos pelo encorajamento e todo
o apoio prestado durante a elaboração desta tese e ao colega Tiago Clara pelo auxílio prestado nesta
fase final do trabalho.
Aos colegas, amigos, que me acompanharam nas inúmeras horas produtivas e não produtivas.
Um especial agradecimento aos meus pais pela paciência e contínuo apoio.
ii
Resumo
Um método potencial de elementos de fronteira (BEM) 3-D foi utilizado para a determinação da matriz
de massa adicionada de um corpo totalmente submerso sujeito a acelerações. A matriz de massa
adicionada é obtida com base nas matrizes de coeficientes de influência que dependem
exclusivamente da malha geométrica que modela o corpo.
Os resultados obtidos serão comparados com valores analíticos e numéricos conhecidos e serão
realizados estudos de convergência com o refinamento de malha, bem como será analisada a
influência da espessura do corpo na qualidade da solução.
À matriz de massa adicionada será ainda aplicado um método de diagonalização e realizado o estudo
da sua capacidade de continuar a fornecer resultados aceitáveis.
Palavras-chave:
Escoamento potencial
Massa adicionada
Método dos painéis
Diagonalização
iii
Abstract
A 3-D potential-based boundary element method is used to compute the added mass matrix of a
totally submerged body subject to accelerations. The added mass matrix is obtained from the
influence coefficient matrices, which are only function of the geometric grid that models the body
surface.
The results obtained are compared against known analytical and numerical values and a grid size
convergence study is conducted. The influence of body thickness in the quality of the solution is
assessed.
A mass diagonalization method is used and the results are compared to the previously obtained
added mass matrix to assess information loss.
Keywords:
Potential flow
Added mass
Panel method
Mass diagonalization
iv
Índice
Agradecimentos ...................................................................................................................... i
Resumo .................................................................................................................................. ii
Palavras-chave:...................................................................................................................... ii
Abstract ................................................................................................................................. iii
Keywords: ............................................................................................................................. iii
Índice .................................................................................................................................... iv
Lista de tabelas ...................................................................................................................... v
1. Introdução ......................................................................................................................... 1
1.1. Introdução ao tema ..................................................................................................... 1
1.2. Trabalhos anteriores ................................................................................................... 1
1.3. Conceito de massa adicionada ................................................................................... 2
1.4. Âmbito do trabalho ...................................................................................................... 3
1.5. Objectivos do presente trabalho .................................................................................. 3
2. Formulação Matemática .................................................................................................... 4
2.1. Massa adicionada ....................................................................................................... 4
2.2. Equação que governa o escoamento .......................................................................... 5
2.2.1. Definição do escoamento ..................................................................................... 5
2.2.2. Equação de Laplace ............................................................................................. 5
2.2.3. Solução geral ....................................................................................................... 7
2.2.4. Equação integral .................................................................................................10
2.2.5 Método dos Painéis ..............................................................................................10
3. Método Numérico .............................................................................................................13
3.1. Geração da malha computacional ..............................................................................13
3.2. Algoritmo do método do painel ...................................................................................14
3.3. Cálculo da matriz de massa adicionada .....................................................................15
3.4. Diagonalização da matriz de massa adicionada.........................................................17
3.5. Cálculo da força hidrodinâmica ..................................................................................17
4. Resultados da Matriz de Massa Adicionada .....................................................................18
4.1 Convergência e verificação .........................................................................................18
4.1.1. Esfera .....................................................................................................................18
4.1.1.1. Esfera – Massa adicionada e respectivo erro ...................................................20
4.1.1.2. Esfera – Visualização das distribuições de forças ............................................23
v
4.1.2. Elipsóide (a=1; b=0,5; c=0,05) ................................................................................28
4.1.2.1. Elipsóide (a=1; b=0,5; c=0,05) – Massa adicionada e respectivo erro ..............29
4.1.2.2. Elipsóide (a=1; b=0,5; c=0,05) – Visualização das distribuições de forças .......32
4.2 Influência da espessura no erro do cálculo da massa adicionada ...............................35
4.3 Diagonalização da matriz de massa adicionada .........................................................41
5. Conclusões ......................................................................................................................43
Referências bibliográficas ....................................................................................................44
Anexo A – Estudo geométrico da malha ...............................................................................45
A.1 – Definição de variáveis ..............................................................................................45
A.2 – Gráficos ...................................................................................................................46
A.2.1 – Esfera ...............................................................................................................46
A.2.2 – Elipsóide (a=1, b=0.5, c=0.05), plano y=0 .........................................................48
Lista de tabelas
Tabela 4-1 - Corpos estudados e suas propriedades geométricas….……………………………………18
Tabela 4-2 – Resultados obtidos para a massa adimensionalizada da esfera e erro relativo quando
esta é sujeita a uma aceleração unitária segundo x………………………………………………………..20
Tabela 4-3 - Resultados obtidos para a massa adimensionalizada da esfera e erro relativo quando
esta é sujeita a uma aceleração unitária segundo y………………………………………………………..20
Tabela 4-4 - Resultados obtidos para a massa adimensionalizada da esfera e erro relativo quando
esta é sujeita a uma aceleração unitária segundo z………………………………………………………..21
Tabela 4-5 - Resultados obtidos para a massa adimensionalizada do elipsóide e erro relativo quando
este é sujeito a uma aceleração unitária segundo x………………………………………………………..29
Tabela 4-6 - Resultados obtidos para a massa adimensionalizada do elipsóide e erro relativo quando
este é sujeito a uma aceleração unitária segundo y………………………………………………………..30
Tabela 4-7 - Resultados obtidos para a massa adimensionalizada do elipsóide e erro relativo quando
este é sujeito a uma aceleração unitária segundo z………………………………………………………..30
Tabela 4-8 - Comparação de parâmetros entre a utilização e não utilização do método de
diagonalização…………………………………………………………………………………………………..41
vi
1
1. Introdução
1.1. Introdução ao tema
Todos os anos, preocupações com alterações climáticas e com o esgotamento dos combustíveis
fosseis colocam mais e mais importância na necessidade da utilização de energias limpas e
renováveis [1,2]. Esta força motriz tem especial impacto na pesquisa e desenvolvimento de
tecnologias de aproveitamento de energia eólica e das marés. Apesar das tecnologias de
aproveitamento de energia das marés produzirem mais energia por unidade de área, é o
aproveitamento de energia eólica que se encontra mais avançada devido a desenvolvimento mais
activo nos últimos 30 anos [2].
À medida que o diâmetro das turbinas eólicas se torna cada vez maior e que os desenvolvimentos de
materiais compósitos permitem a utilização de pás cada vez mais leves, os problemas relacionados
com vibrações tornam-se cada vez mais importantes [1]. Estes mesmos problemas que também se
aplicam a pás propulsoras, nomeadamente de hélices de navios, que devido a variadas razões
operam no seio de fluidos com esteiras não uniformes. Atendendo a que as pás não possuem rigidez
infinita ficam sujeitas a deformações elásticas [3].
Este tipo de problemas provoca uma procura de ferramentas que permitam melhor calcular a
resposta dinâmica e os esforços que se desenvolvem nas pás dos equipamentos de modo a poder
aumentar quer a fiabilidade dos componentes, como reduzir custos de produção e de manutenção
dos mesmos, bem como, no caso de equipamentos produtores de energia poder maximizar a sua
eficiência [3].
1.2. Trabalhos anteriores
Um problema típico de interação fluido-estrutura trata-se da determinação da resposta hidrodinâmica
de uma pá propulsora. Este tipo de problemas pode ser resolvido recorrendo a duas diferentes
famílias de métodos: os métodos de malha de vórtices [VLMs, “Vortex-Lattice Methods”], também
conhecidos por métodos de superfície sustentadora, e por métodos de elementos de fronteira [BEM,
“Boundary Element Methods”] ou método dos painéis.
Os métodos de malha de vórtices foram primeiro utilizados em análises de escoamentos em torno de
propulsores de navios por Kerwin e Lee em 1978 [4]. A primeira utilização dos métodos de elementos
de fronteira para o mesmo tipo de análise foi realizada por Hess e Valarezo em 1985 [5] com base na
velocidade e mais tarde por Lee em 1987 [6] com base no potencial. Apesar de ambos os métodos,
VLM e BEM serem baseados no potencial de velocidades, quando comparados com métodos
2
viscosos, continuam a ser ferramentas de elevada confiança e eficiência para cálculos do
desempenho hidrodinâmico de pás propulsoras [3].
Até recentemente, todos os métodos assumiam que as pás possuem rigidez infinita, sendo que a
resposta estrutural pode ser calculada separadamente usando a distribuição de pressão obtida do
modelo hidrodinâmico. O modelo estrutural mais utilizado é a teoria de viga desenvolvida por Taylor
em 1933 [7].
A teoria de vida demonstrou-se conveniente para estimar os esforços junto da raiz de pás de
geometria convencional, no entanto para pás de geometrias complexas verificou-se que os resultados
não eram aceitáveis. Para fazer face a este problema, a partir da primeira metade da década de 70
começaram a ser usados métodos de elementos finitos [FEMs, “Finite Element Methods”] para as
análises de esforços sobre a estrutura com Genealis em 1970 [8].
Em todos estes trabalhos a pressão do fluido sobre a estrutura era calculada com base na geometria
não deformada das pás. Para considerar as variações na pressão devidas à deformação das pás foi
desenvolvido um método iterativo por Atkinson e Glover em 1988 [9], bem como uma abordagem
acoplada para análise dinâmica dos esforços estruturais por Kuo e Vorus em 1985 [10], em que um
método potencial foi utilizado para determinar as forças hidrodinâmicas bem como a massa
adicionada e amortecimentos associados à deformação de hélices.
1.3. Conceito de massa adicionada
Quando um corpo se move no seio de um fluido, parte do fluido é obrigado a contornar o corpo. Logo,
quando o corpo acelera também o fluido é acelerado o que se traduz na necessidade de exercer mais
força sobre o corpo que a que seria necessária para acelerar o corpo em vácuo. Este excesso de
força é em geral chamado de força adicionada ou força virtual e pode ser contabilizada com a
introdução do conceito de massa adicionada do corpo no seio do fluido.
Considerando a segunda lei de Newton aplicada à translação unidirecional de uma esfera de massa
submerso num fluido é possível escrever:
( ) ( )
onde é a força efectiva sobre o corpo actuado pela aceleração e é a massa adicionada ou
virtual sentida pelo corpo.
3
1.4. Âmbito do trabalho
Para se proceder a uma rigorosa interpretação da interacção fluido-estrutura, torna-se necessário
calcular a matriz de massa adicionada elemento a elemento de modo a permitir tratar de problemas
que envolvam deformação e vibração das pás. Para tal, como se pode verificar em [1,2,3], tem vindo
a ser utilizado um acoplamento de um método de elementos de fronteira com um método de
elementos finitos. O método de elementos de fronteira permite obter a distribuição de pressões e,
consequentemente, a matriz de massa adicionada. Esta, através de um algoritmo de acoplamento, é
adicionada à matriz de massa obtida do método de elementos finitos estrutural (ABAQUS em [1,3]). A
matriz de massa obtida contém a influência do fluido no corpo e pode ser usada em análises
estruturais em modos de vibração.
A mesma metodologia pode ser aplicada para a determinação da matriz de amortecimento
adicionado para a sua adição à matriz de amortecimento estrutural.
Em ambas as situações é possível usar uma técnica de diagonalização para trabalhar com matrizes
com apenas a diagonal principal não nula [3].
1.5. Objectivos do presente trabalho
No presente trabalho será utilizado um método de elementos de fronteira, o método dos painéis, para
a obtenção da distribuição de pressões sobre diversos corpos sujeitos a acelerações. Em seguida
será calculada a matriz de massa adicionada que contabiliza as interações elemento a elemento.
Os objectivos do trabalho podem então ser sintetizados em:
Cálculo de matrizes de massa adicionada para diferentes corpos.
Verificação dos resultados comparando com valores analíticos ou numéricos.
Estudos de convergência com o grau de refinamento da malha de discretização.
Estudos da influência da espessura do corpo nos resultados obtidos.
Aplicação duma técnica de diagonalização à matriz de massa adicionada e comparação no
cálculo das forças com uma matriz não diagonalizada.
4
2. Formulação Matemática
2.1. Massa adicionada
Neste trabalho apenas serão considerados movimentos de translação, nos quais, para um corpo
arbitrário, a massa adicionada não se trata de um escalar como identificada na equação (1.1) mas
sim de uma entidade matricial 3x3, , onde os índices e podem ser x, y ou z:
[ ] [
] ( )
A equação (1.1), apenas para a componente adicionada da massa e força, pode então ser escrita na
forma geral [11]:
( )
onde e representam respectivamente as componentes e dos vectores força e aceleração e
está implícita a soma no índice de 1 a 3.
De (2.1) e (2.2) é fácil perceber que cada componente corresponde à influência na força
devido a uma aceleração .
O conceito de massa adicionada é bastante relevante em estudos de interação fluido-estrutura,
nomeadamente em casos de sistemas mecânicos em vibração num fluido.
A matriz de massa adicionada permite modelar a interação dos elementos da estrutura com as
flutuações de pressão do fluido. Este método tem a vantagem de permitir investigar o comportamento
dinâmico da estrutura sem determinar o campo de escoamentos em torno desta, o que permite
reduzir bastante os graus de liberdade do problema e, consequentemente, o tempo de computação.
Neste trabalho é utilizado o método dos painéis para a determinação das flutuações de pressão do
fluido sobre a estrutura e consequente obtenção da matriz de massa adicionada global de dimensão
3*NP por 3*NP, onde NP representa o número de painéis em que foi dividida a superfície do corpo.
A matriz obtida pode então ser importada para um programa estrutural de elementos finitos, onde
adicionada à matriz de massa da estrutura permite incluir os efeitos do fluido nos cálculos estruturais
[1,2,3].
5
2.2. Equação que governa o escoamento
2.2.1. Definição do escoamento
Considere-se um escoamento em torno de um corpo em movimento no seio de um fluido invíscido e
incompressível em repouso no infinito. Este corpo, ao mover-se, induz no fluido uma perturbação
irrotacional que pode ser representada por um potencial do vector de velocidade.
A velocidade do fluido é dada pela equação (2.3), onde representa o vector velocidades num ponto
e o respectivo potencial.
( )
2.2.2. Equação de Laplace
A equação de continuidade para um escoamento incompressível escreve-se:
( )
A equação que governa o escoamento descrito é a equação diferencial de Laplace para a função
potencial que se obtém da equação da continuidade acima [12].
( ) ( )
Esta equação diferencial linear, em conjunto com as condições de fronteira apropriadas, permite a
solução de escoamentos irrotacionais e incompressíveis em fluido perfeito.
A formulação de condição de fronteira de Neumann implica diretamente que a componente da
velocidade normal à fronteira sólida seja igual à velocidade da fronteira [12].
Para um corpo impermeável submerso num fluido, a componente da velocidade normal às superfícies
sólidas tem de ser nula num referencial fixo com o corpo [12,13,14], isto é:
( )
6
em que é a velocidade relativa da fronteira e a normal unitária apontando para o exterior do
domínio, ver figura 2-1.
A segunda condição de fronteira é a da fronteira a “infinito”, que dita que a velocidade da perturbação
introduzida pelo corpo se reduz a zero com a distância ao mesmo a tender para infinito, equação
(2.8):
( ) ( )
onde, ( ) é o vector de posição em relação ao referencial cartesiano (x,y,z) fixo com o corpo
e é a velocidade relativa entre o corpo e o fluido não perturbado em V.
Figura 2-1 - Identificação da superfície do corpo S e volume
exterior V
7
2.2.3. Solução geral
O método do painel permite a solução de problemas de escoamento potencial. Nestes problemas a
equação da continuidade reduz-se à equação de Laplace como foi acima indicado:
( )
Para obter uma solução a fronteira é descrita por uma linha (em 2-D) ou superfície (em 3-D) divisória
que estabelece a forma do corpo sólido em torno do qual se processa o escoamento. O escoamento
é representado por uma distribuição de singularidades sobre a fronteira [13].
Na base deste método encontra-se o teorema da divergência, apresentado na equação (2.8), que
permite a passagem de integrais de superfície para integrais de volume.
∫
∫
( )
A equação de Laplace (2.5) pode ser resolvida usando a equação (2.8) substituindo o vector q pelo
vector , onde são funções escalares de posição, V é o volume exterior à
superfície S e n a normal apontando para o interior de S [13].
∫ ( )
∫ (
)
( )
Usando a função de Green [3,13],
e tomando , onde r é a distancia de um ponto
P(x,y,z) a outro ponto Q(x’,y’,z’) do domínio de integração e é o potencial do escoamento de
interesse, obtém-se:
∫ (
)
∫
(
) ( )
em que as integrações são realizadas nas variáveis (x’,y’,z’) do ponto Q.
Caso o ponto P seja exterior ao volume V (ver Figura 2-1), então a equação (2.10) toma a forma:
∫ (
)
( )
8
Alternativamente, caso P pertença a V é necessário excluí-lo do domínio de integração, circundando-
o por uma esfera de raio e avaliar o integral sobre a mesma. Nesse caso:
∫ (
)
( )
onde representa a superfície da esfera de raio .
Utilizando coordenadas esféricas e atendendo a que a normal n aponta para o interior da esfera de
raio , isto é, , ⁄ e ( ⁄ ) ( ⁄ ) a equação acima pode ser escrita
como:
∫ (
) ∫ (
)
( )
Avaliando o primeiro integral, fazendo e assumindo que a função e suas derivadas são
contínuas na esfera então a equação (2.12b) pode ser escrita como:
( )
∫ (
)
( )
A equação (2.13a) permite determinar o valor do potencial ( ) para qualquer ponto P interior ao
domínio V em função da influência do potencial ao longo da fronteira S. No caso de P pertencer à
fronteira S, então a integração é apenas efetuada no hemisfério que se encontra em V e (2.12b) toma
a forma:
( )
∫ (
)
( )
Figura 2-2 - Pormenor do
domínio V, ilustrando a
esfera de raio ε em torno
de P
9
É possível também considerar o caso em que o escoamento se processa no interior de S, dando
origem assim a um potencial interno que respeita a equação (2.11):
∫ (
)
( )
Definindo a diferença entre os potenciais exterior e interior, bem como a diferença entre as suas
derivadas, obtém-se uma representação em termos de diferentes tipos de singularidade para a
simulação da superfície do corpo:
(i) Distribuição de dipolos:
( )
(ii) Distribuição de fontes:
( )
Utilizando as expressões (2.15) e (2.16) e atendendo à figura 2-3 é possível reescrever a equação
(2.13a) na forma que será usada ao longo deste trabalho, bem como o seu gradiente:
( )
∫ [ (
)
(
)]
( )
( )
∫ { (
) [
(
)]}
( )
Figura 2-3 – Potencial de
velocidade na vizinhança de uma
parede sólida
10
2.2.4. Equação integral
No caso de ser conhecida a distribuição de singularidades e suas intensidades, então a equação
(2.18) descreve totalmente o campo de velocidades com excepção da fronteira .
Da condição de fronteira (2.6) e tomando o potencial interno como nulo, i.e., bem como
, obtém-se:
( )
e a distribuição de fontes é conhecida. Nesse caso a equação (2.17) para o potencial exterior sobre
escreve-se:
( ) ∫
(
)
∫ (
)
( )
Esta é uma equação integral de Fredholm de 2ª espécie [14], na distribuição de dipolos, que pode ser
resolvida numericamente com o método dos painéis.
2.2.5 Método dos Painéis
No método dos painéis a superfície do corpo é dividida em N painéis e é no centro de cada um
desses painéis que se especificam os pontos de colocação onde a condição de fronteira, através da
equação (2.20), será avaliada.
Para se obter uma solução, a equação acima é avaliada num determinado conjunto de pontos sobre
a superfície do corpo, pontos de colocação, onde para cada um dos pontos é função das
singularidades desconhecidas em todos os pontos. Este método permite reduzir a equação integral
(2.20) a um conjunto de equações algébricas, N equações de N variáveis independentes, sendo N o
número de pontos de colocação sobre a superfície [13].
2.2.5.1. Cálculo dos coeficientes de influência
A equação (2.20), para cada ponto P centrado em cada painel i (ver figura 2-4), assumirá a forma:
∑
∫ { (
) [
(
)]}
( )
11
A equação (2.21) representa o somatório das influências de todos os painéis k. A integração é agora
calculada em cada painel individual e, assumindo singularidades de intensidade unitária, dependerá
apenas da geometria do painel.
Tal como acima mencionado, o problema torma a forma de um sistema de equações algébricas
lineares:
∑
∑
( )
onde e representam os coeficientes de influência do painel k no painel i e são definidos pelas
expressões (2.23a) e (2.23b) respectivamente.
∫
(
) ( )
∫(
) ( )
Figura 2-4 – Discretização da superfície de um disco com identificação dos painéis i e k bem como
do ponto de colocação P do painel i.
12
A equação (2.22) pode ser escrita na forma matricial tomando então a forma [13]:
[
] {
} [
] {
} ( )
Uma vez que é conhecido, ver equação (2.21), e contém a informação sobre a condição de
fronteira de impermeabilidade do corpo sólido, o problema transforma-se na resolução duma equação
matricial completa na forma:
[
] {
} {
} ( )
sendo o vector coluna conhecido, com:
∑
( )
13
3. Método Numérico
Neste capítulo será apresentado o procedimento efetuado para a obtenção da matriz de massa
adicionada.
3.1. Geração da malha computacional
As malhas utilizadas foram geradas recorrendo ao programa PROPANEL, que utilizando uma lei de
coseno, distribui painéis quadriláteros sobre a superfície do corpo de modo a que não haja fugas
entre painéis adjacentes.
A figura 3-1 apresenta um exemplo de uma malha criada pelo programa PROPANEL sobre um
elipsóide de semi-eixos desiguais. A malha possui 200 painéis, e é possível verificar
o efeito do coseno sobre os bordos onde a malha se encontra mais refinada. Na figura 3-1 é o
número de painéis, é o número de painéis segundo a envergadura e é o número de painéis
Figura 3-1 - Discretização de um semi-elipsóide (a=1; b=0,5; c=0,05)
com uma malha 10x10 (NE=10, NC=2x10)
14
segundo a corda. No caso de um elipsóide completo o número de painéis seria o dobro, isto é,
.
Na figura 3-2 é possível verificar diferentes graus de refinamento de malha, respetivamente 400,
1600, 3600 elementos, para uma esfera de raio unitário. É de referir a singularidade da malha nos
polos responsável por erros numéricos acrescidos quando considerado movimento segundo
essa direção.
3.2. Algoritmo do método do painel
O cálculo do potencial de perturbação induzida pela presença do corpo no seio do fluido será
efetuado recorrendo a uma forma discretizada da equação (2.13b) sobre os painéis que simulam a
fronteira do corpo sólido [3]:
∑(
) ( )
Na expressão acima, e , representam respetivamente a influência sentida no painel de
controlo i, devido a um dipolo e uma fonte unitários localizados no painel k.
Figura 3-2 - Da esquerda para a direita, discretização de uma esfera de
raio unitário com as malhas 10x10, 20x20 e 30x30
15
A equação (3.1) pode ser escrita na forma matricial:
[ ]{ } [ ] {
} ( )
Resolvendo em ordem a { } obtém-se
{ } [ ] {
} ( )
com [ ] [ ] [ ]
O mesmo procedimento pode ser usado para as derivadas espaciais [3], isto é:
{
} [ ] {
} ( )
com [ ] [
][ ] [ ], onde [
] e [ ] representam as matrizes de coeficientes da componente
i das velocidades induzidas na direcção , com =1,2,3, por dipolos e fontes unitários,
respetivamente.
Os passos acima mencionados que envolvem o cálculo da inversa da matriz [ ] são os mais
computacionalmente demorados. Neste trabalho utilizou-se um método directo de factorização LU
(sub-rotinas dgefa.f e dgedi.f, Linpack 08/14/78).
3.3. Cálculo da matriz de massa adicionada
Para a determinação da matriz de massa adicionada global, é necessário relacionar painel a painel, a
pressão hidrodinâmica com a área. Para isso comecemos por definir as seguintes igualdades:
[ ]{ } ( )
Definindo o vector de velocidade do centróide do painel k como:
{ } [ ] ( )
onde representam respetivamente as componentes segundo x, y e z do vector de
velocidade do painel k.
Definindo [ ] como o vector das normais do centróide k, obtém-se:
[ ] [
] ( )
16
Atendendo às equações (3.3) e (3.5), o potencial de perturbação no ponto de controlo i pode agora
ser escrito como uma função da matriz de coeficientes de influência, do vector das normais e do
vector de velocidades do painel j, tal como apresentado na equação seguinte:
∑ [ ]{ } ( )
onde NP é o número total de painéis.
Com base na equação acima é possível utilizar uma equação de Bernoulli [3] para obter a pressão
hidrodinâmica sentida no painel i associada à aceleração do fluido de densidade :
( ) ∑ [ ]{ }
( )
Multiplicando a equação (3.9) pelo produto [ ] , onde representa a área do painel i é possível
obter a matriz de massa elementar [ ]:
[ ] [ ]
[ ] ( )
A matriz acima é relativa a um único elemento, possui 3x3 entradas e encontra-se ordenada segundo
as direções x,y,z. Usando esta ordenação é possível proceder à construção [3] da matriz de massa
adicionada global [ ]:
[ ]
[
( )
( )
( )
( )
( ) ( )
]
( )
Cada elemento da matriz (3.11) representa a força hidrodinâmica imposta no grau de liberdade i
por uma aceleração unitária no grau de liberdade j.
A título de exemplo, o elemento da matriz de massa adicionada global, representa a componente
segundo z da força hidrodinâmica exercida no centróide do painel 1 devido a uma aceleração unitária
segundo y do centróide do painel 2.
17
3.4. Diagonalização da matriz de massa adicionada
Um dos objetivos deste trabalho consiste em analisar a possibilidade de utilizar um esquema para
diagonalizar a matriz de massa adicionada e verificar se esta opção cria uma aproximação aceitável
[3].
As entradas da diagonal principal da matriz de massa adicionada diagonalizada [ ] são obtidas
recorrendo à seguinte expressão:
( ) ( ) ( )
∑ ∑ ( ) ( )
∑ ( ) ( )
( )
para e .
O método de diagonalização utilizado garante que a massa total adicionada é conservada
completamente no caso de movimentos de translação [3].
3.5. Cálculo da força hidrodinâmica
Tanto a matriz de massa adicionada global, como a diagonalizada, permitem a obtenção da força
hidrodinâmica através do produto matricial entre a mesma e o vector aceleração a que o corpo se
encontra sujeito, isto é:
[
( )
( )
( )
( )
( ) ( )
]
{
( )
} {
( )
} ( )
onde, por exemplo, e ( ) representam respectivamente a aceleração segundo x do centróide
do painel 3 e a força segundo z no painel NP.
18
4. Resultados da Matriz de Massa Adicionada
Neste capítulo serão apresentados os diversos corpos para os quais foi determinada a matriz de
massa adicionada, bem como, nalguns casos, as forças obtidas e sua distribuição, que servem como
verificação da mesma.
4.1 Convergência e verificação
No presente trabalho foram analisados vários elipsóides, incluindo o caso degenerado da esfera, bem
como uma viga elíptica.
Semi-eixo
Nomenclatura A b C
Esfera 1 1 1
Elipsóide 1 0,5 0,05
Elipsóide 1 1 0,75
Elipsóide 1 1 0,5
Elipsóide 1 1 0,25
Disco 1 1 0,05
Viga elíptica 1* 0,125 0,00625
4.1.1. Esfera
A esfera foi o primeiro corpo a ser analisado devido à fácil obtenção do resultado teórico que permite
validar os resultados numéricos.
A massa adicionada para uma esfera de raio no seio de um fluido de massa específica é dada
pela expressão (4.1) [15]:
( )
Tabela 4-1 – Corpos estudados e suas propriedades geométricas.
* Semi-envergadura
19
Esta expressão pode ser adimensionalizada precisamente pela massa específica do fluido e o cubo
do raio da esfera, tomando a forma:
( )
É com o valor obtido na expressão acima que serão comparados os resultados obtidos pelo método
numérico quando aplicado a uma esfera.
Tal como visto no final do capítulo anterior, a segunda lei de Newton permite utilizar pares
aceleração-força para verificar a matriz de massa adicionada:
( )
sendo e as componentes e dos vectores força e aceleração respectivamente, onde e
podem ser x, y e z e está implícita a soma no índice de 1 a 3.
Figura 4-1 – Quatro malhas criadas pelo programa PROPANEL. De cima para baixo, da
esquerda para a direita, 5x5, 10x10, 20x20, 40x40.
20
4.1.1.1. Esfera – Massa adicionada e respectivo erro
No caso da esfera, foram utilizadas sete malhas diferentes (algumas delas representadas na figura 4-
1), cujo número de elementos, bem como o resultado obtido para a massa adicionada adimensional
são apresentados nas tabelas seguintes para uma aceleração unitária segundo cada uma das
direcções x,y e z.
Definindo a norma de erro ‖ ‖, em que e podem ser x, y e z e que representa o erro relativo
entre a massa adicionada obtida através da equação (4.3) adimensionalizada e o valor analítico
obtido através da expressão (4.2):
‖ ‖ |
|
( )
Nas tabelas 4-2 a 4-4 são apresentados os resultados para a massa adicionada obtidos bem como o
erro definido em (4.4) para as diferentes malhas.
Malha Número de elementos
‖ ‖
5x5 100 2,073037768 7,20E-17 -4,44E-10 1,019737591 10x10 400 2,089768702 -1,46E-16 -6,37E-11 0,220894303 15x15 900 2,092457003 2,28E-16 -2,03E-11 0,092537427 20x20 1600 2,093339874 4,30E-17 -4,22E-12 0,050383424 25x25 2500 2,093733563 8,38E-17 -5,83E-12 0,031586149 30x30 3600 2,093942217 -1,54E-16 -2,96E-12 0,021623675 35x35 4900 2,094065837 4,60E-18 -1,65E-12 0,015721258
Analítico - 2,094395102 0 0 0
Malha Número de elementos
‖ ‖
5x5 100 -1,68E-16 1,904003722 -3,86E-16 9,090518766 10x10 400 4,94E-17 2,047269071 -6,90E-17 2,250102213 15x15 900 -1,25E-16 2,073613417 -1,75E-16 0,992252347 20x20 1600 1,60E-16 2,082760109 3,72E-17 0,555529996 25x25 2500 1,31E-16 2,086971415 3,61E-16 0,354454972 30x30 3600 -4,85E-17 2,089250754 1,43E-16 0,245624550 35x35 4900 -1,63E-16 2,090621518 4,68E-16 0,180175380
Analítico - 0 2,094395102 0 0
Tabela 4-2 – Resultados obtidos para a massa adimensionalizada da
esfera e erro relativo quando esta é sujeita a uma aceleração unitária
segundo x.
Tabela 4-3 – Resultados obtidos para a massa adimensionalizada da
esfera e erro relativo quando esta é sujeita a uma aceleração unitária
segundo y.
21
Malha Número de elementos
‖ ‖
5x5 100 -4,41E-10 -2,12E-16 2,073037779 1,019737097 10x10 400 -6,35E-11 1,39E-16 2,089768714 0,220893766 15x15 900 -2,02E-11 -2,47E-16 2,092457014 0,092536881 20x20 1600 -4,21E-12 6,23E-17 2,093339880 0,050383153 25x25 2500 -5,82E-12 1,78E-16 2,093733579 0,031585421 30x30 3600 -2,95E-12 -2,46E-16 2,093942230 0,021623043 35x35 4900 -1,63E-12 1,16E-16 2,094065847 0,015720792
Analítico - 0 0 2,094395102 0
Dos resultados das tabelas 4-2 a 4-4 é possível verificar que, tal como esperado, o erro diminui
gradualmente com o refinamento da malha. Comparando o caso da aceleração segundo x e da
segundo z, é possível verificar que os valores da massa adicionada adimensionalizada e do
respectivo erro são bastante semelhantes, o que está de acordo com o esperado devido à simetria
total da esfera. No entanto, segundo y o mesmo não se verifica pois é nesse eixo que se encontram
as singularidades do fecho da malha, dando origem a erros manifestamente superiores em
acelerações nessa direcção.
Tabela 4-4 – Resultados obtidos para a massa adimensionalizada da
esfera e erro relativo quando esta é sujeita a uma aceleração unitária
segundo z.
22
1E-2
1E-1
1E0
1E1
1E0 1E1 1E2
||e_xx|| x 100
||e_yy|| x 100
||e_zz|| x 100
No gráfico 4-1 encontram-se representados, em escalas logarítmicas, os erros relativos segundo as
três direcções em função dum parâmetro dependente da malha, a raiz quadrada do número de
elementos.
Este gráfico permite visualizar que o método utilizado é de segunda ordem, isto é, o aumento de uma
ordem de grandeza no parâmetro da malha (raiz quadrada do número de painéis NP) provoca uma
diminuição de duas ordens de grandeza no erro relativo, matematicamente:
‖ ‖ (√ )
( )
A linha sólida no gráfico 4-1 corresponde exatamente a (√ )
.
Gráfico 4-1 – Erros relativos da massa adicionada segundo x quando o corpo é acelerado
segundo x, da massa adicionada segundo y quando o corpo é acelerado segundo y e da
massa adicionada segundo z quando o corpo é acelerado segundo z em função da raiz
quadrada do numero de elementos de malha com ambos os eixos em escala logarítmica.
23
4.1.1.2. Esfera – Visualização das distribuições de forças
Nesta secção serão apresentados vários gráficos onde é possível identificar a distribuição
das forças na superfície da esfera devidas a acelerações unitárias em cada uma das
direcções x, y e z. Malha 30x30, 3600 elementos.
Figura 4-2 – Duas vistas complementares da distribuição da força
segundo x devido a uma aceleração segundo x numa esfera.
Figura 4-3 – Duas vistas complementares da distribuição da força
segundo y devido a uma aceleração segundo x numa esfera.
24
Figura 4-4 – Duas vistas complementares da distribuição da força
segundo z devido a uma aceleração segundo x numa esfera.
Figura 4-5 – Duas vistas complementares da distribuição da força
segundo x devido a uma aceleração segundo y numa esfera.
25
Figura 4-7 – Duas vistas complementares da distribuição da força
segundo z devido a uma aceleração segundo y numa esfera.
Figura 4-6 – Duas vistas complementares da distribuição da força
segundo y devido a uma aceleração segundo y numa esfera.
26
Figura 4-8 – Duas vistas complementares da distribuição da força
segundo x devido a uma aceleração segundo z numa esfera.
Figura 4-9 – Duas vistas complementares da distribuição da força
segundo y devido a uma aceleração segundo z numa esfera.
27
As figuras 4-2 a 4-10 ilustram 18 representações duma esfera discretizada por uma malha 30x30 e as
forças hidrodinâmicas devidas a acelerações do corpo, sendo possível em todas as figuras identificar
a influência das simetrias.
Nas figuras 4-2, 4-6 e 4-10 verifica-se que as maiores intensidades concentram-se junto às calotes
esféricas segundo a direção da aceleração, facto este que se justifica por ser nestes locais que a
perturbação do escoamento torna valores mais elevados. Nestas figuras é também fácil identificar um
certo desvio na concentricidade dos níveis, que se deve à existência da singularidade da malha que
intercepta a esfera nos pontos , responsável também pelos erros numéricos acentuados
quando a aceleração tem essa direcção.
Nas restantes figuras, observa-se como as simetrias provocam o anulamento da resultante das forças
segundo direcções que não a da aceleração. A título de exemplo, na figura 4-9, que representa a
distribuição da força segundo y quando a esfera sofre uma aceleração unitária segundo z, é possível
observar que no tanto no hemisfério como no hemisfério existe uma distribuição
antissimétrica em relação ao plano bissector do eixo z, antissimetria esta que provoca a anulação da
resultante da força em questão.
Figura 4-10 – Duas vistas complementares da distribuição da força
segundo z devido a uma aceleração segundo z numa esfera.
28
4.1.2. Elipsóide (a=1; b=0,5; c=0,05)
O segundo corpo a ser estudado foi um elipsóide de semi-eixos a=1; b=0,5 e c=0,05. Este corpo foi
escolhido por se aproximar da forma geral de uma pá propulsora e assim permitir um cálculo de
potencial mais prático que o corpo anterior, a esfera.
O método para o cálculo da massa adicionada para um elipsóide se eixos arbitrários, a, b e c, foi
obtida de [15]:
( )
onde é a densidade do fluido, são os semi-eixos do elipsóide e é função do eixo que se
encontra alinhado com o movimento:
( )
( )
( )
sendo , e os integrais de Carlson:
∫
( ) ⁄ ( ) ⁄ ( ) ⁄ ( )
∫
( ) ⁄ ( ) ⁄ ( ) ⁄ ( )
∫
( ) ⁄ ( ) ⁄ ( ) ⁄ ( )
que foram resolvidos recorrendo a um programa desenvolvido em [16].
A expressão (4.6) pode ser adimensionalizada utilizando a massa específica do fluido e o cubo do
semi-eixo a, tomando então a forma:
( )
A equação (4.3) serve novamente de base para a obtenção dos resultados a comparar com a solução
numérica obtida de (4.13).
29
4.1.2.1. Elipsóide (a=1; b=0,5; c=0,05) – Massa adicionada e respectivo erro
Tal como no caso da esfera, para o elipsóide também foram utilizadas sete malhas diferentes
(estando algumas representadas na figura 4-11), cujo número de elementos, bem como o resultado
obtido para a massa adicionada adimensional são apresentados nas tabelas seguintes para uma
aceleração unitária segundo cada uma das direcções x,y e z.
Nas tabelas 4-5 a 4-7 são apresentados, para a diferentes malhas, os resultados para a massa
adicionada obtidos bem como o erro, que novamente é definido pela expressão (4.4) mas usando as
variáveis definidas nesta secção para o elipsóide.
Malha Número de elementos
‖ ‖
5x5 100 0,008746044 3,95E-19 -8,70E-11 4,660580221 10x10 400 0,009098586 -4,69E-19 -3,28E-11 0,817567520 15x15 900 0,009144210 5,19E-19 -1,46E-11 0,320226053 20x20 1600 0,009157585 -5,72E-20 -3,73E-12 0,174418116 25x25 2500 0,009163379 2,55E-19 -5,74E-12 0,111264885 30x30 3600 0,009166470 3,69E-19 -3,09E-12 0,077569344 35x35 4900 0,009168333 -4,67E-20 -1,55E-12 0,057260446
Analítico - 0,009173586 0 0 0
Figura 4-11 – Quatro malhas criadas pelo programa PROPANEL. De cima para baixo, da
esquerda para a direita, 5x5, 10x10, 20x20, 40x40.
Tabela 4-5 – Resultados obtidos para a massa adimensionalizada do
elipsóide e erro relativo quando este é sujeito a uma aceleração unitária
segundo x.
30
Malha Número de elementos
‖ ‖
5x5 100 5,99E-19 0,002807324 -2,79E-17 10,84923952 10x10 400 -4,55E-19 0,003063144 8,89E-17 2,725297206 15x15 900 -4,85E-19 0,003111005 -5,81E-17 1,205407668 20x20 1600 -4,53E-20 0,003127663 -1,47E-16 0,676387140 25x25 2500 -2,54E-19 0,003135349 6,50E-17 0,432319190 30x30 3600 -3,64E-19 0,003139517 -2,21E-17 0,299957404 35x35 4900 -6,48E-19 0,003142028 2,15E-16 0,220223884
Analítico - 0 0,003148962 0 0
Malha Número de elementos
‖ ‖
5x5 100 -1,05E-10 6,17E-19 0,944266322 11,14893595 10x10 400 -3,56E-11 -8,30E-19 0,876937135 3,223663912 15x15 900 -1,38E-11 2,30E-18 0,861537771 1,411015381 20x20 1600 -3,13E-12 -4,45E-18 0,856272074 0,791193762 25x25 2500 -4,39E-12 8,91E-19 0,853877372 0,509315045 30x30 3600 -2,21E-12 1,29E-18 0,852578554 0,356432125 35x35 4900 -1,06E-12 2,14E-19 0,851791297 0,263764636
Analítico - 0 0 0,849550483 0
Das tabelas 4-5 a 4-7 é possível verificar que tal como no caso da esfera os erros segundo as três
direcções diminuem com o refinamento da malha. No entanto já não existem as semelhanças de
valores entre os valores para o erro segundo as acelerações nas direcções x e z, isto deve-se ao
facto de o elipsóide apresentar uma espessura muito inferior à da esfera, o que para movimentos
paralelos ao semi-eixo c (ou seja, segundo z) provoca grandes variações do gradiente do potencial de
velocidades em torno do contorno do corpo.
Tabela 4-6 – Resultados obtidos para a massa adimensionalizada do
elipsóide e erro relativo quando este é sujeito a uma aceleração unitária
segundo y.
Tabela 4-7 – Resultados obtidos para a massa adimensionalizada do
elipsóide e erro relativo quando este é sujeito a uma aceleração unitária
segundo z.
31
1E-2
1E-1
1E0
1E1
1E0 1E1 1E2
||e_xx|| x 100
||e_yy|| x 100
||e_zz|| x 100
No gráfico 4-2 encontram-se representados, em escalas logarítmicas, os erros relativos segundo as
três direcções em função dum parâmetro dependente da malha, a raiz quadrada do número de
elementos.
Tal como referido na página anterior, relativamente às tabelas, no gráfico também é visível o aumento
de erro segundo a direcção z devido à inferior espessura do corpo.
Também para o caso do elipsóide, o gráfico permite visualizar que o método utilizado é de segunda
ordem, isto é, o aumento de uma ordem de grandeza no parâmetro da malha (raiz quadrada do
número de painéis NP) provoca uma diminuição de duas ordens de grandeza no erro relativo.
A linha sólida no gráfico 4-2 corresponde exatamente a (√ )
.
Gráfico 4-2 – Erros relativos da massa adicionada segundo x quando o corpo é acelerado
segundo x, da massa adicionada segundo y quando o corpo é acelerado segundo y e da
massa adicionada segundo z quando o corpo é acelerado segundo z em função da raiz
quadrada do numero de elementos de malha com ambos os eixos em escala logarítmica.
32
4.1.2.2. Elipsóide (a=1; b=0,5; c=0,05) – Visualização das distribuições de forças
Nesta secção serão apresentados vários gráficos onde é possível identificar a distribuição das forças
na superfície do elipsóide devidas a acelerações unitárias em cada uma das direcções x, y e z. Malha
35x35, 4900 elementos.
Figura 4-12 – Duas vistas complementares da distribuição da força
segundo x devido a uma aceleração segundo x no elipsóide (a=1; b=0,5;
c=0,05).
Figura 4-13 – Duas vistas complementares da distribuição da força
segundo y devido a uma aceleração segundo x no elipsóide (a=1; b=0,5;
c=0,05).
33
Figura 4-14 – Duas vistas complementares da distribuição da força
segundo z devido a uma aceleração segundo x no elipsóide (a=1; b=0,5;
c=0,05).
Figura 4-15 – Duas vistas complementares da distribuição da força
segundo x devido a uma aceleração segundo y no elipsóide (a=1; b=0,5;
c=0,05).
Figura 4-16 – Duas vistas complementares da distribuição da força
segundo y devido a uma aceleração segundo y no elipsóide (a=1; b=0,5;
c=0,05).
34
Figura 4-17 – Duas vistas complementares da distribuição da força
segundo z devido a uma aceleração segundo y no elipsóide (a=1; b=0,5;
c=0,05).
Figura 4-18 – Duas vistas complementares da distribuição da força
segundo x devido a uma aceleração segundo z no elipsóide (a=1; b=0,5;
c=0,05).
Figura 4-19 – Duas vistas complementares da distribuição da força
segundo y devido a uma aceleração segundo z no elipsóide (a=1; b=0,5;
c=0,05).
35
As figuras 4-12 a 4-20 representam as distribuições de forças sobre o elipsóide quando este é
acelerado segundo as direcções x,y,z. Da observação das mesmas é fácil correlacioná-las com os
valores obtidos para as massas adicionadas apresentadas nas tabelas 4-5 a 4-7.
Tal como no caso da esfera, é possível observar a influência das singularidades da malha quando a
aceleração é segundo y (figura 4-16).
4.2 Influência da espessura no erro do cálculo da massa adicionada
Nesta secção serão apresentados resultados para determinar a influência da espessura na qualidade
dos resultados de massa adicionada. Para tal irá se começar por estudar a esfera, espessura c=1, e
gradualmente fazer decrescer a espessura até se obter um disco de espessura c=0,05. Na figura 4-21
é possível observar os cinco corpos que foram estudados.
A massa adicionada e o respectivo erro para a esfera foram calculados como descrito anteriormente
neste capítulo. No caso dos restantes corpos utilizou-se a metodologia apresentada para o elipsóide
mas como semi-eixos a=1, b=1 e c variável.
Figura 4-20 – Duas vistas complementares da distribuição da força
segundo z devido a uma aceleração segundo z no elipsóide (a=1; b=0,5;
c=0,05).
Figura 4-21 – Corpos utilizados para estudar a influência da espessura.
Da esquerda para a direita: c=1; c= 0,75; c=0,5; c=0,25 e c=0,05.
36
1E-2
1E-1
1E0
1E1
1E0 1E1 1E2
||e_xx|| x 100
||e_yy|| x 100
||e_zz|| x 100
1E-2
1E-1
1E0
1E1
1E0 1E1 1E2
||e_xx|| x 100
||e_yy|| x 100
||e_zz|| x 100
Gráfico 4-3 – Erros relativos da massa adicionada para a esfera (c=1) em função da raiz
quadrada do número de elementos de malha com ambos os eixos em escala logarítmica.
Gráfico 4-4 – Erros relativos da massa adicionada para o segundo corpo da figura 4 -21
(c=0,75) em função da raiz quadrada do número de elementos de malha com ambos os eixos
em escala logarítmica.
37
1E-3
1E-2
1E-1
1E0
1E1
1E0 1E1 1E2
||e_xx|| x 100
||e_yy|| x 100
||e_zz|| x 100
1E-2
1E-1
1E0
1E1
1E0 1E1 1E2
||e_xx|| x 100
||e_yy|| x 100
||e_zz|| x 100
Gráfico 4-5 – Erros relativos da massa adicionada para o terceiro corpo da figura 4 -21 (c=0,5)
em função da raiz quadrada do número de elementos de malha com ambos os eixos em escala
logarítmica.
Gráfico 4-6 – Erros relativos da massa adicionada para o quar to corpo da figura 4-21 (c=0,25)
em função da raiz quadrada do número de elementos de malha com ambos os eixos em escala
logarítmica.
38
Nos gráficos 4-3 a 4-7 a linha sólida é exactamente a mesma que foi definida em (4.5).
Da observação dos cinco gráficos anteriores é verificar que:
1. o erro na massa adicionada segundo x, quando o corpo é acelerado segundo x, aumenta
com a diminuição da espessura,
2. o erro na massa adicionada segundo y, quando o corpo é acelerado segundo y, é
praticamente insensível à diminuição da espessura do corpo, aumenta apenas muito
ligeiramente,
3. o erro na massa adicionada segundo z, quando o corpo é acelerado segundo z, começa por
se reduzir com a diminuição da espessura até ao corpo intermédio (c=0,5) mas, com a
continuação da redução da espessura nos quarto e quinto corpos, acaba por aumentar.
Durante a fase em que o erro diminui, a convergência também se desvia do comportamento
de segunda ordem observado em todos os outros casos.
As duas primeiras observações estão de acordo com o esperado, pois quanto menos suave for a
curvatura maiores serão os gradientes e mais difíceis serão de captar e segundo y existem as
singularidades da malha que induzem erros de magnitude superior que os da variação de espessura.
No entanto, a terceira observação não se encontra de acordo com o esperado, aliás, chega mesmo a
ir contra o que seria intuitivo pois é quando o corpo se move segundo z que os gradientes de
1E-2
1E-1
1E0
1E1
1E0 1E1 1E2
||e_xx|| x 100
||e_yy|| x 100
||e_zz|| x 100
Gráfico 4-7 – Erros relativos da massa adicionada para o quinto e último corpo da figura 4 -21
(c=0,05) em função da raiz quadrada do número de elementos de malha com ambos os eixos
em escala logarítmica.
39
velocidade são mais elevados. O fluido tem de contornar o corpo quase tendo de inverter totalmente
o sentido num percurso extremamente curto (c=0,05).
Para além daquilo que foi exposto no parágrafo anterior, o comportamento da convergência do erro
da massa adicionada segundo z quando o corpo é acelerado nessa mesma direcção também não era
esperado e levou a uma cuidada observação das malhas utilizadas numa tentativa para justificar tal
comportamento. O estudo geométrico das malhas (ver anexo A) revelou que a regra de coseno
utilizada para a geração das malhas é ideal para o caso da esfera pois garante que os painéis sejam
todos semelhantes, mas o mesmo não se verifica para os outros corpos que começam a ter painéis
cada vez mais díspares.
Para a gama de malhas utilizadas, o erro relativo definido na expressão (4.4) da massa adicionada
segundo z quando o corpo é acelerado segundo z, tende para zero por valores negativos para o caso
dos dois primeiros corpos (c=1 e c=0,75) e tende para zero por valores positivos para os dois últimos
corpos estudados (c=0,25 e c=0,05).
-0,02
0
0,02
0,04
0,06
0,08
0,1
0,12
0,14
0,16
10 20 30 40 50 60 70 80
c=1
c=0,75
c=0,50
c=0,25
c=0,05
Gráfico 4-8 – Erro definido pela equação (4.14) para cada um dos cinco corpos em função de
√
40
No caso do corpo intermédio (c=0,5) verificou-se que o erro tende para zero por valores negativos
com o refinamento da malha até 10x10 (400 elementos) e tende para zero por valores positivos para
refinamentos de malha acima de 12x12 (576 elementos). Este facto implica que para um determinado
número de elementos no intervalo ]400,576[ o erro seja nulo.
A justificação acima pode ser constatada nos gráficos 4-8 e 4-9 do erro (4.14) de ⁄ em
função de √ .
( )
Este comportamento dos erros deve-se a erros de integração que se tornam tanto mais importantes
quanto menos uniformes forem os painéis entre si, caso que se verifica quando o corpo se afasta da
esfera.
-0,0005
-0,0003
-0,0001
0,0001
0,0003
0,0005
10 20 30 40 50 60 70 80
c=1
c=0,75
c=0,50
c=0,25
c=0,05
Gráfico 4-9 – Erro definido pela equação (4.14) para cada um dos cinco corpos em função de
√ Ampliação do eixo das ordenadas.
41
4.3 Diagonalização da matriz de massa adicionada
Nesta secção será apresentada uma tabela (Tabela 4-8) onde será possível comparar os resultados
obtidos para a massa adicionada de uma viga elíptica (ver figura 4-22) no caso de se utilizar e de não
se utilizar o método de diagonalização descrito na secção 3.4..
Definindo as variáveis:
(
)
( )
| | |∑∑
∑∑
| ( )
em que (4.15) representa a diferença entre as massas adicionadas segundo x usando e não usando
o método de diagonalização e (4.16) que apresenta a diferença entre a soma de todos os elementos
das matrizes de massa adicionada (com e sem método da diagonalização).
Atendendo às equações (4.15) e (4.16) e aos valores apresentados na tabela 4-8, é fácil de verificar
que para movimentos de translação o método de diagonalização não introduz qualquer tipo de erro
nos resultados obtidos da matriz de massa adicionada diagonalizada quando comparada com
resultados equivalentes da matriz de massa completa.
Malha Número de elementos
(
)
|
|
5x5 100 2,18E-04 2,18E-04 1,00E-18 2,91E-16 10x10 400 2,35E-04 2,35E-04 -3,01E-18 5,55E-17
Tabela 4-8 – Comparação de parâmetros entre a utilização e não
utilização do método de diagonalização.
42
Figura 4-22 – Viga elíptica de semi-envergadura a=1 e semi-eixos da
elipse da secção b=0,125 e c=0,00625.
43
5. Conclusões
Com este trabalho demonstrou-se que o método dos painéis oferece um caminho expedito para o
cálculo de matrizes de massa adicionadas de corpos no seio de fluidos.
Foram obtidos resultados que se aproximam bem aos valores analíticos ou numéricos usados como
verificação e, tal como esperado, observou-se convergência com o grau de refinamento da malha,
sendo a convergência correspondente a um esquema de segunda ordem.
Verificou-se que as singularidades criadas pelo gerador de malhas segundo o eixo y induzem erros
para acelerações segundo y de uma ordem de grandeza superior às outras direcções.
Verificou-se, nomeadamente para o caso de elipsóides, que os erros, segundo as três direcções,
aumentam com a diminuição da espessura do corpo, sendo a direcção z (segundo a espessura) a
mais afectada, seguidamente a direcção x e por ultimo a direcção y, onde a influência da espessura é
quase desprezável devido ao erro induzido pelo referido no parágrafo acima.
De acordo com um dos objectivos propostos, também se provou que o método de diagonalização
utilizado fornece resultados exactos quando se trata de movimentos de translação. O método de
diagonalização é conveniente pois, não só permite uma redução do custo computacional e dos
requerimentos de memória, como transforma a matriz de massa adicionada que geralmente é cheia e
assimétrica numa matriz com apenas a diagonal principal.
A matriz de massa adicionada diagonalizada encontra-se desacoplada para cada grau de liberdade e
pode então ser facilmente adicionada à matriz de massa da estrutura, podendo ser interpretada como
uma distribuição da força inercial do fluido segundo cada direcção nos graus de liberdade que
contribuem para essa direcção.
44
Referências bibliográficas
[1] Young, Y. Motley, M. e Yeung, R., “Three-dimensional numerical modeling of the transient fluid-
structural interaction response of tidal turbines”, J. Offshore Mech. Arctic Eng., Vol. 132. February
2010.
[2] Young, Y. Motley, M. e Yeung, R., “Hydroelastic response of wind or tidal turbines”, 28th
International Conference on Ocean, Offshore and Artic Engineering, Honolulu, Hawaii, USA. 2009.
[3] Young, Y., “Time-dependent hidroelastic analysis of cavitating propulsors”, J. Fluids Struct., 23,
pp. 269-295. 2007.
[4] Kerwin, L., Lee, C.-S., Prediction of steady and unsteady marine propeller performance by
numerical lifting-surface theory. Transactions Society of Naval Architects and Marine Engineers 86,
218-253. 1976.
[5] Hess,J.L., Valarezo, W., Calculation of steady flow about propellers by means of a surface panel
method. In: 23rd
Aerospace Science Meeting. AIAA, Reno, Nevada, pp. 1-8. 1985.
[6] Lee, J.-T., A potencial based panel method for the analysis of marine propellers in steady flow.
Ph.D. Thesis, Department of Ocean Engineering, Massachusetts Institute of Technology, USA.
[7] Taylor, D., The speed and power of ships. Technical report, U.S. Government Printing Office,
Washington, DC. 1933.
[8] Genealis, P., Elastic strength of propellers – an analysis by matrix methods. Technical Report
3397, David Taylor ResearchCenter. 1970.
[9] Atkinson,P., Glover, J., Propeller hydroelastic effects. Transactions Society of Naval Architects
and Marine Engineers 21. 1988.
[10] Kuo, J., Vorus, W., Propeller blade dynamic stress. In: Tenth Ship Technology and Research
(STAR) Symposium. Norfolk, VA, pp. 39-69. 1985.
[11] Munk, M., Fluid Mechanics Pt. 2 in Aerodynamics Theory Vol. 1, edited W. F. Durand, Peter
Smith, 1976.
[12] Mason, B., Applied Computational Aerodynamics.(4. Incompressible Potencial Flow Using Panel
Methods) http://www.dept.aoe.vt.edu/~mason/Mason_f/CAtxtChap4.pdf
[13] Katz, J. e Plotkin, A., Low-Speed Aerodynamics. Singapore. McGraw Hill, 1991. ISBN 0-07-
050446-6.
[14] Brederode, V., Fundamentos de Aerodinâmica Incompressível. Lisboa, Edição de Autor, 1997.
[15] Kennard, E., Irrotational flow of frictionless fluids: Mostly of invariable density, DTMB Report
2299, February 1967.
[16] Oliveira, O., Escoamento Potencial em torno dum elipsóide de eixos desiguais: Solução
Analítica. 1998.
45
Anexo A – Estudo geométrico da malha
Neste anexo será apresentado o estudo geométrico da malha referido na secção 4.2.
A.1 – Definição de variáveis
De modo a encontrar irregularidades nas malhas geradas que pudessem ser a origem do inesperado
comportamento encontrado para o erro segundo z optou-se por definir algumas variáveis que são
função da geometria da malha e observar o seu comportamento com o refinar da malha.
Variáveis para as abcissas:
( )
( )
onde representa o número de pontos da malha na semi-secção em estudo e percorre os pontos
de 1 até .
Variáveis de distância e distância cumulativa:
√( ) ( )
( )
∑ ( )
∑
( )
Razões de espaçamentos:
(
)
( )
(
)
( )
(
)
( )
Na secção seguinte serão apresentados gráficos para as variáveis definidas por (A.3) e (A5-7) em
função de (A.1) e (A.2).
46
A.2 – Gráficos
A.2.1 – Esfera
Devido às simetrias presentes na esfera, os gráficos para o plano y=0 são idênticos aos do plano z=0
pelo que apenas se apresenta os primeiros.
Gráfico A-1 - em função de para a secção com y=0 e z≥0 de uma esfera.
Gráfico A-2 - (
)
em função de para a secção com y=0 e z≥0 de uma esfera.
0
0,1
0,2
0,3
0,4
0,5
0,6
0,7
0,8
0,9
1
0 0,2 0,4 0,6 0,8 1
Esfera 5x5
Esfera 10x10
Esfera 15x15
Esfera 20x20
Esfera 25x25
Esfera 30x30
Esfera 35x35
-2,50E+01
-2,00E+01
-1,50E+01
-1,00E+01
-5,00E+00
0,00E+00
5,00E+00
1,00E+01
1,50E+01
2,00E+01
2,50E+01
0 0,2 0,4 0,6 0,8 1
Esfera 5x5
Esfera 10x10
Esfera 15x15
Esfera 20x20
Esfera 25x25
Esfera 30x30
Esfera 35x35
47
Gráfico A-3 - (
)
em função de para a secção com y=0 e z≥0 de uma esfera.
Gráfico A-4 - (
)
em função de para a secção com y=0 e z≥0 de uma esfera.
0,00E+00
5,00E-01
1,00E+00
1,50E+00
2,00E+00
2,50E+00
3,00E+00
3,50E+00
0 0,2 0,4 0,6 0,8 1
Esfera 5x5
Esfera 10x10
Esfera 15x15
Esfera 20x20
Esfera 25x25
Esfera 30x30
Esfera 35x35
-4,00E+00
-3,00E+00
-2,00E+00
-1,00E+00
0,00E+00
1,00E+00
2,00E+00
3,00E+00
4,00E+00
0 0,2 0,4 0,6 0,8 1
Esfera 5x5
Esfera 10x10
Esfera 15x15
Esfera 20x20
Esfera 25x25
Esfera 30x30
Esfera 35x35
48
A.2.2 – Elipsóide (a=1, b=0.5, c=0.05), plano y=0
Gráfico A-5 - em função de para a secção com y=0 e z≥0 de um elipsóide de semi-eixos a=1,
b=0.5 e c=0.05.
Gráfico A-6 - (
)
em função de para a secção com y=0 e z≥0 de um elipsóide de semi-eixos
a=1, b=0.5 e c=0.05.
0
0,1
0,2
0,3
0,4
0,5
0,6
0,7
0,8
0,9
1
0 0,2 0,4 0,6 0,8 1
5x5
10x10
15x15
20x20
25x25
30x30
35x35
-2,50E+00
-2,00E+00
-1,50E+00
-1,00E+00
-5,00E-01
0,00E+00
5,00E-01
1,00E+00
1,50E+00
2,00E+00
2,50E+00
0 0,2 0,4 0,6 0,8 1
5x5
10x10
15x15
20x20
25x25
30x30
35x35
49
Gráfico A-7 - (
)
em função de para a secção com y=0 e z≥0 de um elipsóide de semi-eixos
a=1, b=0.5 e c=0.05.
Gráfico A-8 - (
)
em função de para a secção com y=0 e z≥0 de um elipsóide de semi-eixos
a=1, b=0.5 e c=0.05.
0
0,2
0,4
0,6
0,8
1
1,2
0 0,2 0,4 0,6 0,8 1
5x5
10x10
15x15
20x20
25x25
30x30
35x35
-1
-0,8
-0,6
-0,4
-0,2
0
0,2
0,4
0,6
0,8
1
0 0,2 0,4 0,6 0,8 1
5x5
10x10
15x15
20x20
25x25
30x30
35x35