Trabalho 2 - Exercício 9 - Mauro Bettoni

Prévia do material em texto

UNIversidade federal do Paraná
MAURO BETTONI JUNIOR
métodos dos elementos finitos
CURITIBA
2021
1. Apresentação da forma fraca para o domínio do problema contendo as condições de contorno.
2. Implementação do problema por meio do pacote FEniCS. Rode o FEniCS e obtenha a evolução temporal da solução até o regime permanente no(s) ponto(s) indicado(s) no problema para intervalos incrementais ∆t, ∆t 2 , ∆t 4 e ∆t 8 e para a família de malhas de elementos lineares: de 2, 4, 8 e 16 elementos. Nos problemas parabólicos, utilize α = 0 e α = 1 2 , e, nos hiperbólicos, utilize α = 1 2 e γ = 1 3 e α = 1 2 e γ = 1 2 .
# -*- coding: utf-8 -*-
"""Untitled0.ipynb
Automatically generated by Colaboratory.
Original file is located at
    https://colab.research.google.com/drive/1iWr6K7qcoYM0Na55KPNfOZWTyrRYH8LF
"""
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 - Trabalho 2
"""
import matplotlib.pyplot as plt
import math
import numpy as np
from dolfin import *
def vibracao(E,k,dt,alpha,gamma) :
    # 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 ponto flutuante
        return on_boundary and abs(x[0]) < tol
    
    # Definição do valor no contorno esquerdo
    u_L = Constant(0)
    
    # definição da malha para problema unidimensional
    mesh = IntervalMesh(E, 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", k)
    
    # Imposição das condições de contorno de Dirichlet à esquerda
    bcs = [DirichletBC(V, u_L, LeftDirichletBoundary)]
    
    n = int(10e-4/dt)
    
    t = []
    u_h = []
    t.append(0)
    
    #Condições iniciais e u e u_vt no instante anterior
    u_s = interpolate(Constant(0), V)
    u_vt_s = interpolate(Constant(0), V)
    u_vtt_s = interpolate(Constant(0), V)
    u_vtt_s1 = interpolate(Constant(0), V)
    
    #Valor de u em x = 0.7
    u_h.append(u_s(0.7))
    
    #Valores do problema 9
    w = Constant(1750e9) # N/m >> Substituto da letra k da mola
    f0 = Constant(3500e3) # N/m²
    r_1 = Constant(0.025) # Raio Menor
    r_2 = Constant(0.05) # 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)
    f1 = Expression('(f0*A_2)/0.4',degree=1,tol=1e-14, f0=f0,A_2=A_2) # N/m
    Eal = Constant(69e9) # GPa =GN/m²
    Eaço = Constant(207e9) # GPa =GN/m²
    L = Constant(0.7)
    
    
    rho_aço = Constant(7980) # kg/m³
    rho_al = Constant(2700) # kg/m³
    EA = Expression('x[0]<=0.4+tol ? Eal*A_2 : x[0]<=0.6+tol ? Eal*A_1 : Eaço*A_1',degree=12,tol=1e-14, Eal=Eal, Eaço=Eaço, A_1=A_1, A_2=A_2)
    rhoA = Expression('x[0]<=0.4+tol ? rho_al*A_2 : x[0]<=0.6+tol ? rho_al*A_1 : rho_aço*A_1',degree=12,tol=1e-14, rho_al=rho_al, rho_aço=rho_aço, A_1=A_1, A_2=A_2)
    
    #Coeficientes
    a = EA
    c2 = rhoA
    f=Expression('x[0]<=0.4+tol ? f0*A_2/0.4:0',degree=2,tol=1e-14,f0=f0,A_2=A_2)
    a1 = Constant((1-alpha)*dt)
    a2 = Constant(alpha*dt)
    a3 = Constant(2/gamma/dt**2)
    a4 = Constant(2/gamma/dt)
    a5 = Constant(1/gamma - 1)
    
    # Definição da forma fraca
    u = TrialFunction(V) # função de aproximação
    v = TestFunction(V) # função peso
    B = a*inner(grad(u),grad(v))*dx + c2*a3*u*v*dx + w*u*v*ds # Operador bilinear
    L = f*v*dx + a3*c2*u_s*v*dx + a4*c2*u_vt_s*v*dx + a5*c2*u_vtt_s*v*dx# operador linear
    
    u = Function(V)
    
    # Loop para sequência temporal
    for i in range(0,n,1):
        # Resolve o problema
        solve( B == L, u, bcs)
        t.append(t[i]+dt)
        u_h.append(u(0.7))
        # Guarda u, u_vt e u_vtt no instante s+1
        u_vtt_s1.assign(a3*(u-u_s) - a4*u_vt_s - a5*u_vtt_s) # u_vtt no instante s+1
        u_vt_s.assign(u_vt_s + a1*u_vtt_s + a2*u_vtt_s1) # u_vt no instante s+1
        u_vtt_s.assign(u_vtt_s1)
        u_s.assign(u)
    print(t)
    print(u_h)
    return u, t, u_h, dt, n, alpha, gamma, E, k
u, t, u_h, dt, n, alpha, gamma, E, k = vibracao(56,2,1e-7,1/2,1/3 ) # tempo(E,k,dt,alpha,gamma)
#plot(u,label='MEF: t={0:.0e} s'.format(t[n]),color ='k')
plt.plot(t,u_h,label='MEF',color = 'k')
# plt.plot(t,u_h,label='FEM',color = 'k')
print(u(0.01))
plt.legend()
plt.title('$\\alpha=${0[0]:,.2f}, $\gamma=${0[1]:,.2f}, $\Delta \, t =${0[2]}, E={0[3]}, k={0[4]}'.format((alpha,gamma,dt,E,k)))
plt.xlabel('x')
plt.ylabel('u(x,{0:.0e})'.format(t[n]))
#plt.xlabel('t')
#plt.ylabel('$u(1,t)$')
3. Apresentação dos resultados do item anterior em gráficos. 
Vibracao(7,2,1e-7,1/2,1/3 )  7 Elementos, dt, gamma = 1/3
vibracao(14,2,1e-7,1/2,1/3 ) 14 Elementos, dt, gamma = 1/3
vibracao(21,2,1e-7,1/2,1/3 ) 21 Elementos, dt, gamma = 1/3
vibracao(28,2,1e-7,1/2,1/3 ) 28 Elementos, dt, gamma = 1/3
vibracao(56,2,1e-7,1/2,1/3 ) 56 Elementos, dt, gamma = 1/3
vibracao(7,2,1e-7/2,1/2,1/3 ) 7 Elementos, dt/2, gamma = 1/3
vibracao(14,2,1e-7/2,1/2,1/3 ) 14 Elementos, dt/2, gamma = 1/3
vibracao(21,2,1e-7/2,1/2,1/3 ) 21 Elementos, dt/2, gamma = 1/3
vibracao(28,2,1e-7/2,1/2,1/3 ) 28 Elementos, dt/2, gamma = 1/3
vibracao(56,2,1e-7/2,1/2,1/3 ) 56 Elementos, dt/2, gamma = 1/3
vibracao(7,2,1e-7/4,1/2,1/3 ) 7 Elementos, dt/4, gamma = 1/3
vibracao(14,2,1e-7/4,1/2,1/3 ) 14 Elementos, dt/4, gamma = 1/3
vibracao(21,2,1e-7/4,1/2,1/3 ) 21 Elementos, dt/4, gamma = 1/3
vibracao(28,2,1e-7/4,1/2,1/3 ) 28 Elementos, dt/4, gamma = 1/3
vibracao(56,2,1e-7/4,1/2,1/3 ) 56 Elementos, dt/4, gamma = 1/3
vibracao(7,2,1e-7/8,1/2,1/3 ) 7 Elementos, dt/8, gamma = 1/3
vibracao(14,2,1e-7/8,1/2,1/3 ) 14 Elementos, dt/8, gamma = 1/3
vibracao(21,2,1e-7/8,1/2,1/3 ) 21 Elementos, dt/8, gamma = 1/3
vibracao(28,2,1e-7/8,1/2,1/3 ) 28 Elementos, dt/8, gamma = 1/3
vibracao(56,2,1e-7/8,1/2,1/3 ) 56 Elementos, dt/8, gamma = 1/3
vibracao(7,2,1e-5,1/2,1/2 ) 7 Elementos, dt, gamma = 1/2
vibracao(14,2,1e-5,1/2,1/2 ) 14 Elementos, dt, gamma = 1/2
vibracao(21,2,1e-5,1/2,1/2 ) 21 Elementos, dt, gamma = 1/2
vibracao(28,2,1e-5,1/2,1/2 ) 28 Elementos, dt, gamma = 1/2
vibracao(56,2,1e-5,1/2,1/2 ) 56 Elementos, dt, gamma = 1/2
vibracao(7,2,1e-5/2,1/2,1/2 ) 7 Elementos, dt/2, gamma = 1/2
vibracao(14,2,1e-5/2,1/2,1/2 ) 14 Elementos, dt/2, gamma = 1/2
vibracao(21,2,1e-5/2,1/2,1/2 ) 21 Elementos, dt/2, gamma = 1/2
vibracao(28,2,1e-5/2,1/2,1/2 ) 28 Elementos, dt/2, gamma = 1/2
vibracao(56,2,1e-5/2,1/2,1/2 ) 56 Elementos, dt/2, gamma = 1/2
vibracao(7,2,1e-5/4,1/2,1/2 ) 7 Elementos, dt/4, gamma = 1/2
vibracao(14,2,1e-5/4,1/2,1/2 ) 14 Elementos, dt/4, gamma = 1/2
vibracao(21,2,1e-5/4,1/2,1/2 ) 21 Elementos, dt/4, gamma = 1/2
vibracao(28,2,1e-5/4,1/2,1/2 ) 28 Elementos, dt/4, gamma = 1/2
vibracao(56,2,1e-5/4,1/2,1/2 ) 56 Elementos, dt/4, gamma = 1/2
vibracao(7,2,1e-5/8,1/2,1/2 ) 7 Elementos, dt/8, gamma = 1/2
vibracao(14,2,1e-5/8,1/2,1/2 ) 14 Elementos, dt/8, gamma = 1/2
vibracao(21,2,1e-5/8,1/2,1/2 ) 21 Elementos, dt/8, gamma = 1/2
vibracao(28,2,1e-5/8,1/2,1/2 ) 28 Elementos, dt/8, gamma = 1/2
vibracao(56,2,1e-5/8,1/2,1/2 ) 56 Elementos, dt/8, gamma = 1/2
4. Comentário dos resultados obtidos.
Conforme vamos aumentando o número de elementos o nível de detalhes dos gráficos apresentados vai diminuindo de forma expressiva, como vemos a seguir:
7 Elementos, DeltaT = 1e-07, gamma = 1/3 _________56 Elementos, DeltaT = 1e-07/4, gamma = 1/3
Quanto mais diminuímos o Δt é necessário muito mais tempo de processamento computacional, utilizando os mesmos parâmetros e equipamentos. Enquantopara um Δt de 0.125e-07 chegou a levar 6 minutos de processamento para gerar a solução e o gráfico, as soluções do Δt a 1e-05 estavam levando cerca de 1 segundo para serem geradas.
Vale a pena notar também que por conta da vibração desta barra em específico, é necessário utilizar sempre de 7 em 7 elementos para que o padrão fosse observado. Caso fosse um número não múltiplo de 7 o gráfico ficaria incoerente, como segue um exemplo com 18 elementos:
Confidential C