Universidade Federal do Espı́rito Santo
Departamento de Informática
Trabalho Computacional de Algoritmos
Numéricos - 15/2
Parte 1: Método de Diferenças Finitas
Data da entrega/apresentação: Enviar o trabalho por email no dia 20/11.
Marcar uma hora para a apresentação do trabalho nos dias 23/11 ou 24/11.
Considerações gerais: O trabalho poder ser feito em grupo de no máximo 3
alunos. Todos os trabalhos devem ser enviados por email (código com todos os arquivos utilizados e o relatório em .pdf compactados. Não inclua os arquivos executáveis
na pasta). O assunto do email deve ”trab1 civil <nome>”, onde ”<nome>”contém
os nomes dos componentes do grupo.
1
Descrição do Problema
Neste trabalho vamos estudar o processo de discretização, pelo método das diferenças
finitas, do modelo que representa o estado estacionário do escoamento de um fluido
ideal em torno de um obstáculo. A aplicação que estamos querendo resolver pode
ser vista como um modelo numérico simplificado do escoamento de um rio raso sobre
um objeto. Neste modelo as hipóteses básicas são: (1) não há viscosidade, ou seja,
não há dissipação de energia devido a atritos internos entre as partı́culas do fluido,
nem das partı́culas com o meio; (2) não há troca de calor entre as partı́culas, nem
das partı́culas com o meio; (3) a densidade é constante em todos os pontos, ou
seja, o fluido é incompressı́vel. Além disso, estamos assumindo que o problema é
bidimensional.
Assumindo um campo de velocidade (u, v), a equação para as linhas de corrente,
definidas por Ψ(x, y), para o modelo mostrado na Fig. 1 é dada por
−
∂ 2Ψ ∂ 2Ψ
−
=0
∂x2
∂y 2
para (x, y) ∈ (0, L) × (0, W )
(1)
satisfazendo as condições de contorno
Ψ
∂Ψ
∂x
Ψ
Ψ
= u0 y
para x = 0 (u = u0 , v = 0)
= 0
para x = L (v = 0)
= 0
= u0 W
para y = 0 ou (x, y) no obstáculo (u = v = 0)
para y = W (u = u0 , v = 0)
onde (u0 , 0) é a velocidade, rio acima, no lado esquerdo do domı́nio e o obstáculo é
definido pelo ponto (XP , YP ). As linhas de corrente são as linhas imaginárias tangentes em cada ponto ao vetor velocidade (u, v). No escoamento estacionário (regime
permanente), as diferentes partı́culas que compõem o fluido têm sempre a mesma
velocidade ao passar em um ponto fixo no espaço. Grandezas como densidade,
velocidade e pressão dependem exclusivamente do ponto do espaço considerado e
independem do tempo. Conseqüentemente, as linhas de corrente irão coincidir com
as trajetórias das partı́culas do fluido. Quando as linhas de corrente se aproximam,
significa que a velocidade do fluido aumentou.
Figura 1: Escoamento em torno de um obstáculo.
2
2.1
Etapas do trabalho
Validação do código com uma solução conhecida
Considere a equação de Poisson
∂ 2u ∂ 2u
+
= f (x, y)
∂x2 ∂y 2
para (x, y) ∈ (0, L) × (0, W )
(2)
satisfazendo as condições de contorno
u
∂u
∂x
u
u
= 25y(y − 1)
x=0
= 0
x=L
= 0
= 0
y=0
y=W
Nesse exemplo a solução exata é conhecida, u(x, y) = (x − 5)2 y(y − 1) com L = 5
and W = 1. Calcule f (x, y) em (2), usando o valor exato da solução dada. São
parâmetros de entrada no seu código: f (x, y), L, W e o número de divisões utilizadas
(N × M ) ou o tamanho das malhas (hx e hy ) nas direções x e y. Além disso, tem que
ser fornecido o parâmetro w do método SOR que será utilizado, um número máximo
de iterações SOR (maxSOR) e uma tolerância para o método SOR (tolSOR).
Nos testes computacionais: (1) faça um estudo do parâmetro ω ∈ (0, 2) para o
método SOR para diferentes malhas: 20 × 4, 40 × 8, 80 × 16 e 160 × 32. Para cada
malha, mostre em um gráfico ou em uma tabela a variação de ω com o número de
iterações necessários para a convergência (um dos parâmetros deve ser w = 1) para
a escolha de um w “ótimo”; (2) Faça também um estudo da exatidão da solução
calculando a norma do erro dada por: erro = max|uexato
−ui | para i = 1, 2, ·, N ∗M .
i
Mostre em uma tabela a malha e o erro máximo calculado. Escolha uma malha e
mostre os gráficos das soluções exata e aproximada obtidas, considerando a solução
pelo método SOR com o ω “ótimo” encontrado.
2.2
Escoamento em torno de um obstáculo
Nos testes computacionais do escoamento em torno de um obstáculo, considere o
obstáculo indentificado pelo ponto (xP , yP ). Considere os seguintes dados para o
problema:
1. u0 = 1, L = 500, W = 100, N = 50 e M = 20, (XP , YP ) = (400, 70),
maxSOR = 1000 e tolSOR = 0.01. Escolha o parâmetro w, considerando
os experimentos para o problema com solução conhecida. Compare com o
número de iterações para o caso de w = 1;
2. Resolva também considerando o obstáculo localizado em (XP , YP ) = (300, 50).
Traçe o gráfico dos contornos das linhas de corrente para cada caso e comente os
resultados. Utilize o octave, o gnuplot ou outro aplicativo para fazer os gráficos.
Para isso, salve os resultados em um arquivo de saı́da.
3
Estrutura do relatório
O relatório dever ser escrito observando as normas do padrão ABNT. A divisão do
relatório deve ser de acordo com as seguintes seções:
• Introdução: onde o grupo dever apresentar a estrutura do trabalho e os
objetivos;
• Método das Diferenças Finitas: onde o grupo descrever as equações do
modelo matemático e a discretização por diferenças finitas utilizada.
• Experimentos Numéricos onde serão apresentados os testes computacionais, com as entradas para os programas, as tabelas e gráficos das respectivas
saı́das geradas pelas soluções e as discusões dos resultados obtidos.
• Conclusão: onde é apresentado um resumo do que foi feito no trabalho juntamente com as conclusões.
• Apêndice: onde serão apresentados a estutura do código e partes significativas do código comentado.
Parte 2: EDO usando Octave/Matlab
Utilizar o octave/matlab para resolver equações diferenciais ordinárias. Coloque no
seu código todos o comentários e resultados, de forma a facilitar a apresentação.
Exemplos de como deve ser feita a estrutura do trabalho encontram-se na página do
curso.
Experimentos Computacionais:
1. Implemente o método de Runge-Kutta de ordem 4 no octave/matlab para
um problema com solução conhecida e verifique a ordem de convergência do
método.
2. Considere o problema de valor inicial
dy
2
2
+ 0.6y = 10e−(x−2) /[2(0.075) ]
dx
y(0) = 0.5
(3)
(4)
Calcule a solução no intervalo [0, 4], com h = 0.1 e h = 0.05, mostrando em
um mesmo gráfico as duas soluções. Observe que a maior discrepância entre
os dois resultados ocorre na região entre 1.8 e 2.0.
3. Utilize o método de Runge-Kutta Adaptativo ode45 para resolver o problema,
mostrando a solução em um gráfico. Mostre também os passos escolhidos pelo
algoritmo.
4. A equação diferencial básica de uma curva elástica para uma viga em balanço
(Figura 2) é dada por
d2 y
EI 2 = −P (L − x)
dx
onde E é o módulo de elasticidade e I é o momento de inércia. Determine
a deflexão da barra usando um método numérico. Os seguintes valores dos
parâmetros se aplicam: E = 30.000 ksi, I = 800 in4 , P = 1 kip, L = 10 pés.
x3
Lx2
Compare seus resultados numéricos com a solução analı́tica, y = − P2EI
= P6EI
,
mostrando as soluções em um mesmo gráfico.
Figura 2: Viga em balanço.