LACTEA - Computation Technology Laboratory of Environmental Engineering

Search

Equilibrio.py

[codesyntax lang=”python”]

# -*- coding: utf-8 -*-
"""
Data: 21/11/2016
@Autor: Prof. Emílio G. F. Mercuri D.Sc.
Departamento de Engenharia Ambiental (DEA)
Disciplina: TEA008 - Mecanica dos Solidos II
Curso: Engenharia Ambiental UFPR

Trabalho: Dinâmica de Vibração de Árvores
Programa que resolve a equação de equilíbrio

Referência: Murphy, K. D., & Rudnicki, M. (2012). 
A physics-based link model for tree vibrations. 
American journal of botany, 99(12), 1918-1929.

"""

print("""
Trabalho: Dinamica de Vibracao de Arvores
Programa que resolve a equacao de equilibrio
""")


# Importação das Bibliotecas e Funções
import numpy as np
from scipy import optimize
from scipy import linalg
from scipy.optimize import fmin
from scipy.optimize import fsolve

# Definição das variáveis
# Versão atual: caso simétrico (1 tronco e 2 ramos)
grav   = 9.81 # m/s2
L      = 6.1  # m
m      = 43.8 # kg
K      = 162.7*10**3 # N.m
h1     = 4.3 # m
h2     = 4.3 # m
theta1 =   np.pi/4 # rad
theta2 = - np.pi/4 # rad
L1     = 3.7 # m
L2     = 3.7 # m
m1     = 14.6 # m
m2     = 14.6 # m
K1     = 8.8*10**3 # N.m
K2     = 8.8*10**3 # N.m


# Definição das funções da Matriz M

def P1(x):
    return m*(L**2)/3 + m1*h1*(0.5*L1*np.cos(x[1])*np.cos(theta1)+h1-0.5*L1*np.sin(x[1])*np.sin(theta1)) + m2*h2*(0.5*L2*np.cos(x[2])*np.cos(theta2)+h2-0.5*L2*np.sin(x[2])*np.sin(theta2))

def P21(x):
    return  m1*h1*0.5*L1*(np.cos(theta1)*np.cos(x[1])-np.sin(theta1)*np.sin(x[1]))

def P22(x):
    return  m2*h2*0.5*L2*(np.cos(theta2)*np.cos(x[2])-np.sin(theta2)*np.sin(x[2]))

def Q11(x):
    return  m1*((0.5*L1)**2)*np.sin(theta1+x[1])*(np.cos(x[1])*np.sin(theta1)+np.sin(x[1])*np.cos(theta1)) + m1*0.5*L1*np.sin(0.5*np.pi-theta1-x[1])*(h1+0.5*L1*np.cos(x[1])*np.cos(theta1)-0.5*L1*np.sin(x[1])*np.sin(theta1)) + (1/12)*m1*(L1)**2   

def Q12(x):
    return  m2*((0.5*L2)**2)*np.sin(theta2+x[2])*(np.cos(x[2])*np.sin(theta2)+np.sin(x[2])*np.cos(theta2)) + m2*0.5*L2*np.sin(0.5*np.pi-theta2-x[2])*(h2+0.5*L2*np.cos(x[2])*np.cos(theta2)-0.5*L2*np.sin(x[2])*np.sin(theta2)) + (1/12)*m2*(L2)**2 

def Q21(x):
    return  m1*((0.5*L1)**2)*np.sin(theta1+x[1])*(np.cos(x[1])*np.sin(theta1)+np.sin(x[1])*np.cos(theta1)) + m1*((0.5*L1)**2)*np.sin(0.5*np.pi-theta1-x[1])*(np.cos(x[1])*np.cos(theta1)-np.sin(x[1])*np.sin(theta1)) + (1/12)*m1*(L1)**2   

def Q22(x):
    return  m2*((0.5*L2)**2)*np.sin(theta2+x[2])*(np.cos(x[2])*np.sin(theta2)+np.sin(x[2])*np.cos(theta2)) + m2*((0.5*L2)**2)*np.sin(0.5*np.pi-theta2-x[2])*(np.cos(x[2])*np.cos(theta2)-np.sin(x[2])*np.sin(theta2)) + (1/12)*m2*(L2)**2   
    
# Montagem da matrix M
def M(x):
    return np.matrix([[P1(x) , P21(x), P22(x)],
                      [Q11(x), Q21(x), 0],
                      [Q12(x), 0     , Q22(x)]])    
                      
x = np.zeros(3)
#print(M(x))


# Definição das componentes do vetor g(x)

def g1(x):
    return K*x[0] - K1*x[1] - K2*x[2] - m*grav*0.5*L*np.sin(x[0]) - (h1*m1*grav*np.sin(x[0])+h2*m2*grav*np.sin(x[0]))

def g2(x):
    return K1*x[1] -m1*grav*0.5*L1*np.sin(theta1+x[1])*np.cos(x[0]) - m1*grav*0.5*L1*np.sin(0.5*np.pi-theta1-x[1])*np.sin(x[0])

def g3(x):
    return K2*x[2] -m2*grav*0.5*L2*np.sin(theta2+x[2])*np.cos(x[0]) - m2*grav*0.5*L2*np.sin(0.5*np.pi-theta2-x[2])*np.sin(x[0])

def g(x):
    return np.array([[g1(x)],[g2(x)],[g3(x)]]) 
    
#print(g(x))

# Cálculo da inversa M-1(x)
#print(linalg.inv(M(x)))   

# Cálculo do vetor f (como um vetor coluna 1x3)
def f(x):
    return linalg.inv(M(x)).dot(g(x))
    
#print(f(x))

# Cálculo do vetor f (como um vetor linha)
def func(x):
    out = [linalg.inv(M(x)).dot(g(x))[0][0]]
    out.append(linalg.inv(M(x)).dot(g(x))[1][0])
    out.append(linalg.inv(M(x)).dot(g(x))[2][0])
    return out
 
#print(func(x))
 
# Chute inicial para resolução do sistema de equações não linear
x0 = [0.0, 0.7, 0.1]
# Resolução da eq. do equilíbrio com a função fsolve (scipy)
solucao = fsolve(func, x0)

# imprime a solução em radianos
print("Solução do equilibrio em radianos:")
print(solucao)


#Convertendo a solução de radianos para graus:
thetachapeu0 = solucao[0]*180/np.pi
thetachapeu1 = solucao[1]*180/np.pi
thetachapeu2 = solucao[2]*180/np.pi


# imprime a solução em radianos
print("")
print("Solução do equilibrio em graus:")
print("")
print("Theta chapeu tronco: %.4f graus" % thetachapeu0)
print("Theta chapeu ramo 1: %.4f graus" % thetachapeu1)
print("Theta chapeu ramo 2: %.4f graus" % thetachapeu2)

[/codesyntax]

Back to Top