Trabalho 1 - Exercício 9 - Mauro Bettoni

Prévia do material em texto

Confidential C 
UNIVERSIDADE FEDERAL DO PARANÁ 
 
 
MAURO BETTONI JUNIOR 
 
 
 
 
 
 
 
 
 
 
MÉTODOS DOS ELEMENTOS FINITOS 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
CURITIBA 
2021 
 
 
 
Confidential C 
 
 
 
 
 
Confidential C 
 
1. A forma fraca do problema para um elemento genérico (e). 
 
 
 
2. O modelo local de elementos _nitos para os elementos linear e quadrático 
genéricos (para os problemas que envolvem viga, apresente o modelo só para o 
elemento de viga). 
 
 
 
 
Confidential C 
 
 
 
 
 
Confidential C 
 
 
 
 
 
 
Confidential C 
 
 
 
 
Confidential C 
 
 
 
 
 
 
Confidential C 
3. Os respectivos modelos globais de elementos finitos (linear e quadrático) para 
malha de 2 elementos com as condiicões de contorno impostas (para os problemas 
que envolvem viga, apresente o modelo só para o elemento de viga). 
 
 
 
 
 
 
 
 
 
Confidential C 
4. Implemente o problema por meio do pacote FENICS. Rode o FENICS para a 
família de malhas de elementos lineares 1 (2, 4, 8, 16, 32, 64, 128 e 256 elementos) 
e para a família de malhas de elementos quadráticos (2,4,8,16,32,64,128 e 256 
elementos). Apresente juntos, para comparação, os gráficos da solução exata e das 
soluções linear e quadrática com 8 elementos. 
 
try: 
 import dolfin 
except ImportError: 
 !wget "https://fem-on-colab.github.io/releases/fenics-install.sh" -
O "/tmp/fenics-install.sh" && bash "/tmp/fenics-install.sh" 
 import dolfin 
 
""" 
Exercício 9, Mauro Bettoni Junior 
 
""" 
import matplotlib.pyplot as plt 
import math 
import numpy as np 
from dolfin import * 
 
# definição da malha para problema unidimensional 
mesh = IntervalMesh(256, 0, 0.7) 
 
# Definição do espaço de função na malha utilizando polinômios de Lagra
nge do 1o. grau 
V = FunctionSpace(mesh, "Lagrange",1) 
 
# Condição de Dirichlet à esquerda 
def LeftDirichletBoundary( x, on_boundary): 
 # Função que verifica se ponto x pertence ao LeftDirichletBoundary 
com tolerânica de +/- tol 
 # retorna veradeiro ou falso 
 tol = 1E-14 
 # tolerância necessária ao lidar com representação de número em pon
to flutuante 
 return on_boundary and abs(x[0]) < tol 
 
# Definição do valor no contorno esquerdo 
u_L = Constant(0) 
 
# Imposição das condições de contorno de Dirichlet à esquerda 
bcs = [DirichletBC(V, u_L, LeftDirichletBoundary)] 
 
# Definição da forma fraca 
Eal = Constant(69e9) # Pa = N/m² 
Eaço = Constant(207e9) # Pa = N/m² 
w = Constant(1750e9) # N/m >> Substituto da letra k da mola 
 
 
 
Confidential C 
L = Constant(0.7) # m 
P = Constant(450e3) # N 
f0 = Constant(3500e3) # N/m² 
f = Expression('((f0*A_2)-(2*P))/0.4',degree=1,tol=1e-
14, f0=f0,A_2=A_2,P=P) # N/m 
gamma = Constant(0.7) # da condição de Robin 
r_1 = Constant(0.0125) # Raio Menor 
r_2 = Constant(0.025) # Raio Maior 
pi = Constant(3.14159265358979323846) 
A_1= Expression('pi*r_1*r_1',degree=2, tol=1e-14, pi=pi,r_1=r_1) 
A_2= Expression('pi*r_2*r_2',degree=2, tol=1e-14, pi=pi,r_2=r_2) 
 
u = TrialFunction(V) # função de aproximação 
v = TestFunction(V) # função peso 
a = Expression('x[0]<=0.4+tol ? Eal*A_2 : x[0]<=0.6+tol ? Eal*A_1 : Eaç
o*A_1',degree=2,tol=1e-14, Eal=Eal, Eaço=Eaço, A_1=A_1, A_2=A_2) 
B = a*inner(grad(u), grad(v))*dx - (w-f)*u*v*ds # operador bilinear 
L = f*v*dx - (w-f)*gamma*v*ds # operador linear 
 
# Resolve o problema 
u = Function(V) 
solve( B == L, u, bcs) 
 
# Plota a solução 
plot(u,label='256 elementos - Linear') 
 
plt.legend() 
 
plt.title('Exercício 9') 
 
 
 
 
 
 
 
Confidential C 
 
 
 
 
 
 
 
Confidential C 
 
 
 
 
 
 
 
Confidential C 
5. Para cada família do item anterior, levante o gráfico em escala logarítmica do erro 
nas normas L2 e da energia em função do número de graus de liberdade e estime a 
taxa de convergência para cada família. Compare cada estimativa com o valor 
teórico esperado. 
 
try: 
 import dolfin 
except ImportError: 
 !wget "https://fem-on-colab.github.io/releases/fenics-install.sh" -
O "/tmp/fenics-install.sh" && bash "/tmp/fenics-install.sh" 
 import dolfin 
 
""" 
Exercício 9, Mauro Bettoni Junior 
 
""" 
import matplotlib.pyplot as plt 
import math 
import numpy as np 
from dolfin import * 
 
# Condição de Dirichlet à esquerda 
 
def LeftDirichletBoundary(x1, on_boundary): 
 # Função que verifica se ponto x pertence ao LeftDirichletBoundary 
com tolerânica de +/- tol 
 # retorna veradeiro ou falso 
 tol = 1E-14 
 # tolerância necessária ao lidar com representação de número em pon
to flutuante 
 return on_boundary and abs(x1) < tol 
 
# Definição do valor no contorno esquerdo 
u_L = Constant(0) 
 
# Solução exata 
 
Eal = Constant(69e9) # N/m² 
Eaço = Constant(207e9) # N/m² 
w = Constant(1750e9) # N/m 
L = Constant(0.7) # m 
P = Constant(450e3) # N 
f = Constant(3500e3) # N/m² 
gamma = Constant(1) # da condição de Robin 
r_1 = Constant(0.0125) # Raio Menor 
r_2 = Constant(0.025) # Raio Maior 
pi = Constant(3.14159265358979323846) 
 
 
 
Confidential C 
A_1= Expression('pi*r_1*r_1',degree=2, tol=1e-14, pi=pi,r_1=r_1) 
A_2= Expression('pi*r_2*r_2',degree=2, tol=1e-14, pi=pi,r_2=r_2) 
 
alfa = Expression('(Eal*A)/(w*L)',degree=2, tol=1e-
14,Eal=69e9, A=4.91e-4, w=1750e9,L=0.7) 
beta = Expression('(P)/(f*L*L)',degree=2, tol=1e-
14, P=450e3, f=3500e3,L=0.7) 
 
u_e = Expression('x[0]<=0.4+tol ? -(1/392)*((f*x[0]*(-336*(L*L)*alfa-
384*(L*L)+4116*(L*L)*beta*alfa+343*(x[0]*x[0])*alfa+490*(x[0]*x[0])))/(
alfa*w*(7*alfa+10)*L)) : x[0]<=0.6+tol ? -(2/343)*(f*L*(-
8+147*beta)*(19*L+7*L*alfa-21*x[0]))/(alfa*w*(7*alfa+10)) : -
(2/49)*(f*L*(-8+147*beta)*(L*alfa+L-x[0]))/(alfa*w*(7*alfa+10))', 
 degree=2, tol=1e-
14, alfa=alfa, beta=beta, f=3500e3, w=1750e9,L=0.7) 
 
def Erro(k, faixa_1, faixa_2): 
 # Array para os erros de cada malha 
 E_L2 = [] 
 E_H1 = [] 
 # Loop para sequência de malhas 
 for i in faixa_1: 
 # definição da malha para problema unidimensional 
 mesh = IntervalMesh(2**i, 0, 0.7) 
 
 # Definição do espaço de função na malha utilizando polinômios 
de Lagrange do 1o. grau 
 V = FunctionSpace(mesh, "Lagrange", 1) 
 
 # Imposição das condições de contorno de Dirichlet à esquerda e
 à direita 
 bcs = [DirichletBC(V, u_L, LeftDirichletBoundary)] 
 
 # Definição da forma fraca 
 u = TrialFunction(V) # função de aproximação 
 v = TestFunction(V) # função peso 
 a = Expression('x[0]<=0.4+tol ? Eal*A_2 : x[0]<=0.6+tol ? Eal*A
_1 : Eaço*A_1',degree=2,tol=1e-
14, Eal=Eal, Eaço=Eaço, A_1=A_1, A_2=A_2) 
 B = a*inner(grad(u), grad(v))*dx - w*u*v*ds # operador bilinear 
 L = f*v*dx - w*gamma*v*ds # operador linear 
 
 # Resolve o problema 
 u = Function(V) 
 solve(B == L, u, bcs) 
 
 # Erro na norma L² 
 
 
 
Confidential C 
 E_L2.append(errornorm(u_e, u, 'L2')) 
 # Erro na norma H¹ 
 E_H1.append(errornorm(u_e, u, 'H1')) 
 # Norm = norm(u_e,'H1',mesh) 
 # print('Norma =', Norm) 
 # print('Erro relativo =', E/Norm*100,'% ') 
 # Taxas de convergência 
 x = [1/2**i for i in faixa_1] # x = { 1, 1/2, 1/2³, 1/2⁴, ..., 1/2¹
⁰} 
 p_L2 = [math.log(E_L2[i]/E_L2[i-1])/math.log(x[i]/x[i-
1]) for i in faixa_2] 
 p_H1 = [math.log(E_H1[i]/E_H1[i-1])/math.log(x[i]/x[i-
1]) for i in faixa_2] 
 return x, E_L2, E_H1, p_L2, p_H1 
 
# Definindo os pares ordenado (x,y) e (x,z) para plotagem do gráfico lo
g-log do erro 
x, y, z, p, q = Erro(2, range(0, 11, 1), range(1, 11, 1)) # Norma L² e
 Norma H¹ 
 
print(p) 
print(q) 
 
plt.plot(x, y, label='L²') # O gráfico 
plt.plot(x, z, label='H¹') # O gráfico 
plt.xscale("log") # Escala logarítmica em x 
plt.yscale("log") # Escala logarítmica em y 
 
plt.title('k= 1') 
 
plt.xlabel('h') 
plt.ylabel('erro') 
 
#plt.annotate('local max', xy=(1e-3, 1e-9), xytext=(1e-2, 1e-
9), arrowprops=dict(facecolor='black', shrink=0.05)) 
 
plt.grid(True) 
 
plt.legend() # mostra a legenda 
 
# Definindo o par ordenado (x,y) para plotagem do gráfico de u_exata 
# x = np.linspace( 0, 1, 50 ) 
# y = [a+2.5/math.pi**2*math.sin(2*math.pi*a) for a in x] 
 
# Plota a solução 
 
 
 
 
Confidential C