Modelação Numérica 2016/17

Exercícios de revisão

Tópico 2. Equações diferenciais parciais.

Mestrado Integrado em Engenharia da Energia e do Ambiente
Licenciatura em Meteorologia, Oceanografia e Geofísica
Licenciatura em Engenharia Geoespacial

DEGGE-FCUL
Departamento de Engenharia Geográfica, Geofísica e Energia
Faculdade de Ciências da Universidade de Lisboa

Conteúdos

  1. Diferenças finitas: avançadas, retardadas e centradas.
  2. Métodos implícitos e explícitos.
  3. Métodos de Lax, Upstream, Leapfrog.
  4. Condições fronteira abertas, fechadas e periódicas.
  5. Estabilidade.
  6. Método de relaxamento de Jacobi.
  7. Problemas de condições iniciais e de valores de fronteira.

Materiais de estudo

A referência principal para estudo são os slides das aulas teóricas.

Como material auxiliar, diponibilizam-se também os anotamentos da cadeira desenvolvidos pelo Prof. Fernando Santos. Note-se que apesar da grande sobreposição entre as aulas deste ano e os anotamentos, a sobreposição não é completa. Os anotamentos cobrem algumas matérias que não foram dadas este ano (e que, portanto, não serão avaliadas). Por outro lado, também houve tópicos cobertos nas aulas que não estão nos anotamentos (e que poderão ser avaliados). Notas sobre equações diferenciais.

Sugestões

No exame não haverá acesso a computador. No entanto, para estudar, pode ser vantajoso utilizar o computador para manipular os códigos aqui apresentados, bem como os códigos apresentados nas aulas, para melhor os perceber.

Exercício 1.

Considera uma nuvem de poluente que, na ausência de vento, se vai difundir. Para modelar este fenómeno, utilizamos a equação da difusão:

$$ \frac{\partial \chi}{\partial t} = -D \frac{\partial^2 \chi}{\partial x^2} $$

1. Discretiza a equação acima, utilizando um método explícito. Identifica que tipo de diferenças utilizaste para discretizar cada derivada parcial (retardadas, avançadas ou centradas).

2. Como poderias alterar a discretização acima, de forma a torná-la num esquema implícito?

3. Esquematiza o essencial do código (pseudo-código) que utilizarias para modelar a difusão da nuvem de poluente (não é necessário indicar as linhas de código em detalhe, basta identificar os blocos de código que seriam necessários e para que servem).

4. Repete o exercício 1, adimitindo agora que a difusão ocorre nas três dimensões do espaço:

$$ \frac{\partial \chi}{\partial t} = -D \left[ \frac{\partial^2 \chi}{\partial x^2} + \frac{\partial^2 \chi}{\partial y^2} + \frac{\partial^2 \chi}{\partial z^2} \right] $$

5. Se em vez de ser difundida, a nuvem de poluente for apenas advectada (transportada) pela força do vento, a equação que descreve a sua evolução é a equação da advecção:

$$ \frac{\partial \chi}{\partial t} = -u \frac{\partial \chi}{\partial x} $$

Discretiza a equação acima utilizando os métodos de Lax-Friedrichs, Upstream e Leapfrog. Como é que estes métodos diferem de um simples método "Forward in time, centered in space"?

6. Como farias a implementação de condições fronteira periódicas? E abertas? Em que condições utilizarias cada um deste tipo de fronteiras?

Exercício 2.

Considera o campo potencial electromagnético gerado por 3 cargas eléctricas pontuais (caso 2D):

$$ \frac{\partial^2 V}{\partial x^2} + \frac{\partial^2 V}{\partial y^2} = -\frac{\rho}{\epsilon} $$

O código abaixo, bem como o seu output, modela este campo potencial utilizando o método iterativo de Jacobi

In [2]:
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from math import pi as pi
%matplotlib inline

plt.rcParams['figure.figsize'] = 10, 6

#%% Parâmetros
dx=.01                   
dy=.01                   
dx2dy2=dx**2/dy**2
a=2.*(1+dx2dy2)

x=np.arange(0.,1.,dx); nx=len(x)      
y=np.arange(0.,1.,dy); ny=len(y)      

V=np.zeros([nx,ny])         
f=np.zeros([nx,ny])         

c1=1e-9
c2=-1e-9
c3=5e-10
eps=8.88541878172e-12

f[40,50]= -c1/eps
f[60,50]= -c2/eps
f[50,15]= -c3/eps

it=0            

Vold=np.zeros([nx,ny]) 

V[0,:]=V[1,:]       
V[nx-1,:]=V[nx-2,:]
V[:,0]=V[:,1]
V[:,ny-1]=0

#%% Plot 3D 
fig = plt.figure()
ax = plt.axes(projection='3d')

yy, xx = np.meshgrid(x, y)         
surf=ax.plot_surface(xx,yy, V, rstride=5, cstride=5, cmap=cm.jet, linewidth=0)                 
fig.colorbar(surf, shrink=0.5)

ax.set_xlabel(u'Posição x (m)')
ax.set_ylabel(u'Posição y (m)')
ax.set_zlabel(u'Potencial Eléctrico')
Out[2]:
<matplotlib.text.Text at 0x10e5e1790>
In [3]:
#%% Plot 3D 
fig = plt.figure()
ax = plt.axes(projection='3d')

surf=ax.plot_surface(xx,yy, f, rstride=5, cstride=5, cmap=cm.jet, linewidth=0)                 
fig.colorbar(surf, shrink=0.5)

ax.set_xlabel(u'Posição x (m)')
ax.set_ylabel(u'Posição y (m)')
ax.set_zlabel(u'Forçamento f')
Out[3]:
<matplotlib.text.Text at 0x111105290>
In [5]:
#%% Método de Jacobi

niter = 1000        

for i in range(niter):
    it += 1
    Vold[:]=V[:]    
    ms=0            
    for ix in range(1,nx-1):
        for iy in range(1,ny-1):
            V[ix,iy] = (V[ix+1,iy] + V[ix-1,iy] + dx2dy2*(V[ix,iy+1] + V[ix,iy-1]) \
                        - dx**2*f[ix,iy]) / a
                        
    V[0,:]=V[1,:]       
    V[nx-1,:]=V[nx-2,:]
    V[:,0]=V[:,1]
    V[:,ny-1]=V[:,ny-2]

#%% Plot 3D
           
fig = plt.figure()
ax = plt.axes(projection='3d')

surf=ax.plot_surface(xx,yy, V, rstride=5, cstride=5, cmap=cm.jet, linewidth=0)                 
fig.colorbar(surf, shrink=0.5)

ax.set_xlabel(u'Posição x (m)')
ax.set_ylabel(u'Posição y (m)')
ax.set_zlabel(u'Potencial Eléctrico')
ax.set_title(u'Iteração = '+str(it))
Out[5]:
<matplotlib.text.Text at 0x10bafed10>

1. Discretiza a equação diferencial do campo potencial electromagético, utilizando diferenças centradas no espaço.

2. Qual a relação entre a discretização que fizeste e a equação implementada no código?

3. Como alterarias o código se quisesses modelar o campo potencial electromagnético gerado apenas pelas cargas positivas?

4. E como farias se quisesses adicionar uma carga negativa em y=0.8 m e x = 0.9 m? De que informação adicional precisarias?

5. Como modificarias o código, se precisasses de maior resolução espacial?

6. Como modificarias o código, se precisasses de modelar um domínio espacial maior?

7. Que tipo de condições fronteira é utilizado no código acima?

8. Como modificarias o código, se o campo electromagnético na fronteira do domínio estivesse fixo em valores nulos?

9. Como podes avaliar se o problema já foi iterado num número de vezes suficientes para a solução estar correcta?

In [ ]: