Vibrações - Multiplos Gdl

PEC00025 githubmmaiarochaPEC00025treemaster Classe16TestP2dipynb githubmmaiarochaPEC00025treemasterClass16TestP2dipynb Universidade Federal do Rio Grande do Sul UFRGS Programa de PósGraduação em Engenharia Civil PPGEC PEC00025 Introdução à Teoria da Vibração Teste P2 20231 Sistemas mdof discretos e contínuos NOME CARTÃO Em 1 Em 2 Questão 1 O sistema estrutural abaixo com 2 gdl representa um pórtico plano dotado de um amortecedor de massa sintonizada Calcular os coeficientes de instabilidadek 1 ek 2 e as respectivas formas modais correspondentes às frequências naturais dadas Importando módulos Python necessários para este notebook esta célula deve ser executada com shiftenter antes de qualquer outra célula P importar numpy como np importar matplotlibpyplot como plt importar pickle como pk importar scipylinalg como sc de MRPy importar def modosdevibração K M Usa scipy para resolver o problema de autovalor padrão w2 Phi sc eig K M Garanta a ordem crescente dos autovalores iw w2 argsort w2 w2 iw Phi Phi iw Autovalores para frequências de vibração wk np sqrt np real w2 fk wk 2 np pi retornar fk wk Phi Solução As matrizes de massa e de sofrimento são K k 1 k 2 k 2 k 2 k 2 M M 1 0 0 M 2 Lembrando agora o problema de autovalores e autovetores K φ eu ω2 eu M φ eu Considerando que as formas modais são normalizadas pela sua coordenada emM 1 k 1 k 2 k 2 k 2 k 2 1 φ eu ω2 eu M 1 0 0 M 2 1 φ eu Isso resulta em um par de equações k 1 k 2 k 2 φ eu ω2 euM 1 k 2 k 2 φ eu ω2 euM 2 φ eu páraeu 1eeu 2 Portanto temse um sistema com 4 equações para 4 incógnitas k 1 k 2 φ 1 φ 2 A dificuldade está na nãolinearidade Somandose as duas equações k 1 ω2 euM 1 M 2 φ eu Isolandosek 2 na primeira equação k 2 ω2 euM 1 k 1 1 φ eu Inicialmente vamos esperar um valor inicial às vezes para se obter as formas modais e então tentar iterar Em 3 Frequência no primeiro modo 089 Hz Frequência no segundo modo 112 Hz Rigidez das colunas inferiores k1 394784 Nm Rigidez do TMD k2 19739 Nm Coordenada do primeiro modo em M2 500 Coordenada do segundo modo em M2 400 Após algumas tentativas percebese que a única maneira de se obter as frequências dadas é reduzindo a massa M2 para 100kg 1 da massa M1 ou aumentandose a massa M1 para 50 ton mesma relação As rigidezes são calculadas usandose a média das duas frequências alvo que é 1Hz Portanto adotase as frequências 089 e 112Hz obtidas acima correspondentes às rigidezes propostas Questão 2 Para a estrutura do problema anterior calcule as máximas amplitudes de deslocamento de cada massa para uma carga dinâmica harmônica aplicada na massa M1 Ft F0sin2πf0t M1 10000 M2 500 w12 2nppi1002 w22 2nppi1002 k1 w12M1 k2 w22M2 KG nparray k1k2 k2 k2 k2 MG nparray M1 0 0 M2 fk wk Phi vibrationmodesKG MG ph1 Phi0Phi00 normalizando pela coordenada da massa M1 ph2 Phi1Phi01 printFrequência no primeiro modo 062f Hzformatfk0 printFrequência no segundo modo 062f Hz formatfk1 printRigidez das colunas inferiores k1 060f Nmformatk1 printRigidez do TMD k2 060f Nm formatk2 printCoordenada do primeiro modo em M2 062fformatph11 printCoordenada do segundo modo em M2 062fformatph21 onde F0 1kN e f0 1Hz Método 1 solução numérica por Duhamel Inicialmente construímos o vetor de cargas NODAIS sendo uma força harmônica na massa M1 e zero na massa M2 In 4 Cálculo dos parâmetros modais usando a matriz Φ original fornecida pelo algoritmo de autovalores do Python com a escala que tiver Abaixo também são calculadas as forças modais In 5 Simulação das forças NODAIS Td 128 N 216 F0 1000 f0 100 t nplinspace0 Td N F1 F0npsin2nppif0t F2 npzeroslikeF1 FG MRPynpvstackF1 F2 TdTd FGplottimefig1 figsize126 axist0 FGTd 2000 2000 zk nparray001 001 Mk npdiagnpdotPhiT npdotMG Phi Kk Mkwk2 Fk MRPynpdotPhiT FG fsFGfs Fkplottimefig2 figsize126 axist0 FGTd 500 500 Cálculo dos deslocamentos por superposição MODAL e retorno aos valores NODAIS In 6 Deslocamento modal de pico do primeiro modo 287mm Deslocamento modal de pico do segundo modo 231mm Deslocamento de pico da massa M1 11mm Deslocamento de pico da massa M2 505mm ak MRPynpdotnpdiag1Mk Fk fsFkfs divide force by modal mass uk aksdofDuhamelfk zk modal space solution uG 1000MRPynpdotPhi uk fsukfs back to nodal solution uG uGextractsegm23 1 byfraction avoid transiente start uGplottime4 figsize126 axist0 uGTd 60 60 f uGfaxis Su fs uGperiodogram pltfigure5 figsize124 pltsemilogyf Su0 b f Su1 g pltaxis0 5 1e4 1e5 pltlegendMassa M1 Massa M1 pltgridTrue ukp ukextractsegm34 1 byfractionmaxaxis1 up uGmaxaxis1 printDeslocamento modal de pico do primeiro modo 041fmmformat1000ukp0 printDeslocamento modal de pico do segundo modo 041fmm format1000ukp1 printDeslocamento de pico da massa M1 041fmmformatup0 printDeslocamento de pico da massa M2 041fmmformatup1 Observase que a frequência da carga está situada bem entre as duas frequências naturais de modo que se evita um forte pico de ressonância Método 2 pela função de admitância As propriedades modais já foram calculadas no método anterior Agora calculase a amplitude das forças modais In 7 Amplitude da força modal no primeiro modo 1961 N Amplitude da força modal no segundo modo 2425 N A amplitude dos deslocamentos modais é calculada usandose a função de ganho Aβ ζ na frequência da excitação f0 1Hz FG nparrayF0 0T amplitudes das forças nodais vetor coluna Fk npdotPhiT FG Fk1 Fk00 Fk2 Fk10 printAmplitude da força modal no primeiro modo 041f NformatFk1 printAmplitude da força modal no segundo modo 041f NformatFk2 In 8 Amplificação dinâmica no primeiro modo 398 Deslocamento modal no primeiro modo 2859 mm Amplificação dinâmica no segundo modo 498 Deslocamento modal no segundo modo 2312mm Observase que os deslocamentos MODAIS coincidem com os valores obtidos por simulação Finalmente os deslocamentos MODAIS são superpostos ignorandose a fase para se calcular os deslocamentos NODAIS Como a informação de fase é perdida é feita uma combinação quadrática das amplitudes que traz alguma imprecisão In 9 Deslocamento de pico da massa M1 79mm Deslocamento de pico da massa M2 359mm Observase que os resultados diferem dos obtidos por simulação onde a fase é levada em conta Questão 3 Para a viga com as restrições de apoio dadas proponha uma forma aproximada para o primeiro modo de vibração φx e calcule a respectiva frequência natural em função do comprimento L da massa por unidade de comprimento μ e da rigidez à flexão EI Aw lambda f npsqrt1 1 ffn22 2ztffn2 fn fk0 zt zk0 A1 Awf0 u1 1000A1Fk1Kk0 fn fk1 zt zk1 A2 Awf0 u2 1000A2Fk2Kk1 printAmplificação dinâmica no primeiro modo 052fformatA1 printDeslocamento modal no primeiro modo 052f mm formatu1 printAmplificação dinâmica no segundo modo 052fformatA2 printDeslocamento modal no segundo modo 052fmmformatu2 uG npsqrtPhi0u12 Phi1u22 Combinação quadrática de printDeslocamento de pico da massa M1 041fmmformatuG0 printDeslocamento de pico da massa M2 041fmmformatuG1 Solução Vamos usar como modo de vibração a função de interpolação para uma rotação do nó B Logo ϕξ ξ3 ξ2 onde ξ xL A escala da função acima não é importante A curvatura é a segunda derivada dessa função ϕξ 6ξ 2L2 Logo a energia cinética de referência é 2Tref 1 0μϕ2ξ Ldξ μL 105 e a energia potencial elástica é 2V 1 0EI ϕξ 2 Ldξ 4EI L3 Portanto pelo quociente de Rayleigh temos fn 1 2π V Tref 1 2π 420EI μL4 1 2π 45270 L 2 EI μ Na tabela abaixo o resultado exato é encontrado como sendo 393 ao invés de 453 com a função aproximada Lembrando que o quociente de Rayleigh sempre dá uma frequência acima do valor correto O gráfico abaixo é só para conferir a forma modal aproximada escolhida In 10 Alternativamente podemos usar a linha elástica que resulta da aplicação de uma carga distribuída sobre uma viga com as condições de contorno dadas ϕξ 2ξ4 5ξ3 3ξ2 phi lambda xi Lxixixi xixi L 1 x nplinspace0 L 1024 phx phix pltfigure6 figsize124 pltplotx phx pltaxis0 1 02 01 pltgridTrue onde ξ xL A escala da função acima é irrelevante A curvatura é a segunda derivada dessa função ϕξ 24ξ2 30ξ 6 L2 Logo a energia cinética de referência é 2Tref 1 0μϕ2ξ Ldξ 19μL 630 e a energia potencial elástica é 2V 1 0EI ϕξ 2 Ldξ 36EI 5L3 Portanto pelo quociente de Rayleigh temos fn 1 2π V Tref 1 2π 22680EI 95μL4 1 2π 39308 L 2 EI μ Observase que essa função de interpolação proposta é quase exatamente o valor apresentado na tabela Abaixo está um gráfico para visualização da forma modal proposta In 11 Questão 4 Para a viga do problema anterior com os valores dados abaixo calcule a amplitude máxima de deslocamento para o peso próprio sendo aplicado de forma súbita a partir do tempot 0 ou seja como uma função passo unitárioht phi lambda xi 2xixixi 5xixixi 3xixi L 1 x nplinspace0 L 1024 phx phix pltfigure6 figsize124 pltplotx phx pltaxis0 1 06 02 pltgridTrue Solução Inicialmente definase os dados numéricos do problema Em 12 Usandose a função de interpolação proposta são calculados os valores modais Primeiro a massa modal M k 1 0μϕ 2 ξ eudξ 19μeu 630 452kg e em seguida a frequência modal c k 3931 eu 2 EEU μ 1881rumds 2993Hpor com as quais calculamos a rigidez modal K k ω2 kM k 1600kNm A amplitude da força modal passo unitário é calculada a partir da forma modal F k 1 0μgϕξ eudξ 3 20μgeu 221kN A amplificação dinâmica para uma carga passo unitária éUM 2 que deve ser aplicado sobre o deslocamento modal estático vocêk den UMvocêk est UM F k K k 2 221 1600 276mm A amplitude máxima da forma modal é calculada como sendoϕmumx 026 Portanto a amplitude máxima de deslocamento é dada por L 6 comprimento das barras m EI 48000000 sofrimento à flexão Nm2 μ 250 massa por unidade de comprimento kgm g 981 gravidade ms2 q μ g carga por unidade de comprimento Nm você máx ϕmumxvocêk den 072mm Em PEC00025 githubmmaiarochaPEC00025treemaster Classe16TestP2cipynb githubmmaiarochaPEC00025treemasterClass16TestP2cipynb Universidade Federal do Rio Grande do Sul UFRGS Programa de PósGraduação em Engenharia Civil PPGEC PEC00025 Introdução à Teoria da Vibração Teste P2 20221 Sistemas mdof discretos e contínuos NOME CARTÃO Instruções 1 Entregar a resolução da prova em arquivo único com no máximo 10Mb até às 24h de hoje 25 de maio de 2022 2 Recomendase verificar atentamente se todas as folhas de resolução foram incluídas no arquivo gerado pois não serão aceitas entregas posteriores 3 Na primeira folha do arquivo deverá constar claramente o NOME e o cartão de MATRÍCULA 4 A consulta ao material de estudo e o uso do computador para cálculos são LIVRES 5 A prova deve ser realizada INDIVIDUALMENTE sem recorrer ao auxílio de colegas ou outras pessoas Caso se verifique o descumprimento desta regra todos os envolvidos na fraude terão uma nota da prova zerada Em 1 Importando módulos Python necessários para este notebook esta célula deve ser executada com shiftenter antes de qualquer outra célula P importar numpy como np importar matplotlibpyplot como plt importar pickle como pk importar scipylinalg como sc de MRPy importar MRPy Em 2 Questão 1 Um cabo com comportamento elástico linear é disposto horizontalmente O cabo tem comprimento totaleu 6m rigidez axialEUM 4000kN e duas massasm 20kgintroduzidos nos terços O cabo tem uma tensão inicialT 0 20kN O amortecimento do sistema éζ 001 razão do crítico A flexibilidade à flexão bem como a massa do cabo são desprezíveis A gravidade da gravidade no local ég 981ms 2 Os dois graus de liberdade considerados são os deslocamentos verticais das duas massasvocê 1 tevocê 2 t Solução Admitindose uma condição de pequenos deslocamentos calcule os modos de vibração e as respectivas frequências naturais de vibração livre do sistema def modosdevibração K M Usa scipy para resolver o problema de autovalor padrão w2 Phi sc eig K M Garanta a ordem crescente dos autovalores iw w2 argsort w2 w2 iw Phi Phi iw Autovalores para frequências de vibração wk np sqrt np real w2 fk wk 2 np pi retornar fk wk Phi In 3 Questão 2 No grau de liberdade u1t do problema anterior é aplicada uma carga transiente F1t dada pela função abaixo com amplitude F0 500N e duração Td 01s A variável τ representa o tempo adimensionalizado por Td L 60 m 20 T0 20000 k 6T0L K1 nparrayk k2k2 k M1 nparraym 0 0 m fk1 wk1 Phi1 vibrationmodesK1 M1 f1 pltfigure1 figsize105 x nparange0 8 2 for k in range2 qk npzeros4 qk11 Phi1k qk npmaxnpabsqk adjust scale for unity amplitude pltsubplot21k1 pltplotx qk pltxlim 00 60 pltylim15 15 pltylabelstrk1 plttext35 06 fk 042fHzformatfk1k fontsize14 pltgridTrue pltxlabelx Desconsiderando a parcela estática da resposta devida ao peso próprio e considerando todos os modos de vibração apresente o deslocamento u1t como uma função do tempo Indique a amplitude e o instante no tempo em que o máximo deslocamento é atingido Solução Método 1 por superposição modal simulando o carregamento In 4 Simulação das forças NODAIS N 8192 Td 01 F0 500 t nplinspace0 20Td N τ tTd fs Nt1 F1 4F0τ τ2 F1t Td 0 F1 MRPynpvstackF1 npzeroslikeF1 fsfs f2 F1plottimefig2 figsize104 axist0 F1Td 200 800 In 5 In 6 Deslocamento de pico da massa 1 é 00426m Deslocamento de pico da massa 2 é 00378m Cálculo das forças MODAIS zk1 nparray001 001 Mk1 npdiagnpdotPhi1T npdotM1 Phi1 Kk1 Mk1wk12 Fk1 MRPynpdotPhi1T F1 fsF1fs f3 Fk1plottimefig3 figsize105 axist0 F1Td 500 500 Cálculo dos deslocamentos por superposição modal Mass division must be a matrix operation ak1 MRPynpdotnpdiag1Mk1 Fk1 fsFk1fs and now solving uk1 ak1sdofDuhamelfk1 zk1 modal space solution uN1 MRPynpdotPhi1 uk1 fsuk1fs back to nodal solution Resultado no domínio do tempo f4 uN1plottime4 figsize105 axist0 uN1Td 008 008 Resultado no domínio da frequência para confirmar picos f5 uN1plotfreq5 figsize105 axisf0 10 0 0001 up1 uN1maxaxis1 printDeslocamento de pico da massa 1 é 064fmformatup10 printDeslocamento de pico da massa 2 é 064fmformatup11 Método 2 por condições iniciais carga como impulso de Dirac In 7 Mode 1 with phase 000rad and amplitude 00527m Mode 2 with phase 000rad and amplitude 00304m In 8 Deslocamento de pico da massa 1 é 00577m Deslocamento de pico da massa 2 é 00579m I0 2F0Td3 impulse is Ftdt v I0m which is converted to a initial veloci u0 nparray0 0T column vector with the initial displacemen v0 nparrayv 0T column vector with the initial velocities qMu npdotnpdotPhi1T M1 u0 qMv npdotnpdotPhi1T M1 v0 thk npzeroslikeMk1 phase angles to be calculated u0k npzeroslikeMk1 modal response amplitude to be calculated for k in range2 If there are initial displacements only thkk nppi2 u0kk qMukMk1knpsinthkk If there are initial velocities only thkk nparctanwk1kqMukqMvk u0kk qMvkMk1knpcosthkkwk1k printMode 0 with phase 152frad and amplitude 274fmformatk1 thk Build the modal responses as harmonic functions with given properties uk MRPyharmonicNX2 NN fsF1fs X0u0k f0fk1 phithk Calculate the NODAL responses superposing all modal responses uN MRPynpdotPhi1 uk fsF1fs f6 uNplottime6 figsize105 axist0 uNTd 008 008 up uNmaxaxis1 printDeslocamento de pico da massa 1 é 064fmformatup0 printDeslocamento de pico da massa 2 é 064fmformatup1 Observase que a hipótese de aproximação como carga impulsiva é muito ruim pois para o segundo modo td Tn 01s 616Hz 0616 e este valor é bem maior que o valor de 25 recomendado como limite para a aproximação Além disso devese considerar que a solução acima não considera o amortecimento Questão 3 Todos os elementos do pórtico elástico linear tem rigidez à flexão EI 65kNm2 e massa por unidade de comprimento μ 20kgm Viga e colunas tem o mesmo comprimento L 4m O amortecimento do sistema é ζ 001 razão do crítico A aceleração da gravidade no local é g 981ms2 Proponha funções adequadas para representar uma geometria deformada que aproxime o primeiro modo de vibração e estime a frequência fundamental de vibração livre através do quociente de Rayleigh Lembre que as energias totais serão computadas somandose a contribuição dos três elementos estruturais Sugerese o uso do software Ftool para o cálculo da energia interna de deformação Solução Fazendose o cálculo da estrutura acima no Ftool temse a seguinte deformada Os deslocamentos nos nós superiores calculados com o Ftool com uma carga estática F 1kN são In 9 Estes valores podem ser interpolados utilizando as funções de interpolação dadas em aula para termos uma expressão analítica para a linha elástica com a qual podese calcular a energia de deformação e a energia cinética de referência In 10 uA 0590 deslocamento da extremidade superior esquerda m uB uA deslocamento da extremidade superior direita m θA 00894 rotação da extremidade superior esquerda rad θB 00887 rotação da extremidade superior direita rad Dados do problema L 4 comprimento das barras m EI 6500 rigidez à flexão Nm2 μ 20 massa por unidade de comprimento kgm F 1000 carga estática arbitrária aplicada N Discretização do comprimento das barras x nplinspace0 L 200 dx L200 Lambda functions para interpolação dos deslocamentos phi phiappendlambda xi 1 3xixi 2xixixi phiappendlambda xi Lxi 2xixi xixixi phiappendlambda xi 3xixi 2xixixi phiappendlambda xi Lxixi xixixi Lambda functions para interpolação das curvaturas phixx phixxappendlambda xi 6 12xiLL phixxappendlambda xi 4 6xiL phixxappendlambda xi 6 12xiLL phixxappendlambda xi 2 6xiL In 11 Esse conjunto de linhas deformadas será utilizado como forma modal para o cálculo da respectiva frequência natural de vibração livre através do quociente de Rayleigh Observe que além dos deslocamentos transversais a cada barra temse também o deslocamento da viga para a direita que representa a maior parte da energia cinética do sistema In 12 A energia cinética de referência é 2258 Jm A massa da viga para a direita representa 617 A energia potencial elástica pode ser calculada pelo trabalho da força externa Deslocamentos interpolados para as três barras coluna da esquerda x de baixo pra cima w1 0phi0xL 0phi1xL uAphi2xL θAphi3xL viga superior x da esquerda para direita w2 0phi0xL θAphi1xL 0phi2xL θBphi3xL coluna direita x de baixo pra cima w3 0phi0xL 0phi1xL uBphi2xL θBphi3xL f7 pltfigure7 figsize55 s 2 escala das deformações pltplotsw1 x b suA x L sw2 b L sw3 x b pltaxisequal pltgridTrue Energia cinética de referência Deslocamento da viga para a direita Tv μLuA22 Deslocamentos transversais das três barras Tr Tv μnptrapzw12 w22 w32 dxdx2 printA energia cinética de referência é 042f JmformatTr printA massa da viga para a direita representa 031fformat100TvTr In 13 A energia potencial elástica é 2950 J E finalmente o cálculo pelo quociente de Rayleigh In 14 A frequência fundamental do pórtico é menor que 0575 Hz O cálculo da energia potencial elástica também pode ser feito pela curvatura In 15 A energia potencial elástica é 2955 J Que resulta muito próximo do valor calculado pelo trabalho das forças externas respeitando portanto a conservação de energia Questão 4 O topo do pórtico é submetido a uma força horizontal estocástica Ft com densidade espectral SFf ilustrada abaixo A força tem média zero e valor rms σF 50N A banda de frequências excitada é definida por r 10 eixo das frequências em hertz Trabalho da força externa V FuA2 printA energia potencial elástica é 041f JformatV fn npsqrtVTr2nppi printA frequência fundamental do pórtico é menor que 053f Hzformatfn Curvaturas w1xx 0phixx0xL 0phixx1xL uAphixx2xL θAphixx3xL w2xx 0phixx0xL θAphixx1xL 0phixx2xL θBphixx3xL w3xx 0phixx0xL 0phixx1xL uBphixx2xL θBphixx3xL V EInptrapzw1xx2 w2xx2 w3xx2 dxdx2 printA energia potencial elástica é 041f JformatV Estime o valor rms e o valor de pico do deslocamento horizontal ut na extremidade esquerda da viga Calcule a correspondente força estática equivalente Solução Inicialmente vamos calcular as propriedades modais no primeiro e único modo de vibração que foi estimado Neste caso a configuração deformada aproxima a forma modal cuja escala vamos manter com o deslocamento horizontal uA na extremidade esquerda da viga Definida essa escala para a forma modal a massa modal iguala a energia cinética de referência sem o fator 12 In 16 Massa modal 452 kg Rigidez modal 5900 Nm Deslocamento estático 0590 m A força definida pelo espectro é aplicada na extremidade esquerda da viga na horizontal onde o deslocamento na forma modal tem amplitude uA Portanto o espectro da força modal é Cálculo das propriedades modais Mk 2Tr wk 2nppifn Kk wkwkMk zk 001 Deslocamento estático a partir da rigidez modal Fk 1000uA força modal uk FkKk deslocamento estático modal ue uAuk deslocamento nodal printMassa modal 051f kgformatMk printRigidez modal 041f NmformatKk printDeslocamento estático 053f mformatue In 17 Amplitude rms pela integral do espectro é 500 N Uma vez definido o espectro da força modal podemos calcular a resposta modal no domínio da frequência O espectro do deslocamento NODAL extremidade esquerda da viga é obtido multiplicando se o espectro do deslocamento MODAL por u2 A da mesma forma como foi feito para a força modal M 4097 discretização do domínio da frequência σF 50 valor rms da força r 10 f nplinspace0 2r M fs 2f1 SF npzeroslikef SFf 1r σF22nplogr1ff 1r SFf 1r 0 SFf r 0 sF2 nptrapzSF f sF npsqrtsF2 printAmplitude rms pela integral do espectro é 041f N formatsF SFk uAuASF espectro da força modal Fkt pltfigure3 figsize84 clearTrue pltsemilogyf SFk pltgridTrue pltaxis0 12r 1e01 1e04 pltxlabelFrequência Hz pltylabelEspectro da força modal N²Hz In 18 Finalmente fazse a análise estatística da resposta em deslocamento a partir do espectro In 19 O fator de pico da resposta em deslocamento é 287 A amplitude rms da resposta no eixo é 123mm O valor do pico da resposta em deslocamento é 353mm Podemos usar a aparência aparentemente da análise do Ftool para calcular a força estática equivalente Em 20 Força estática equivalente é 5982N Em Hf2 lambda fi 1 1 fifn22 2zkfifn2 Kk2 SU uAuAHf2fSFk pltfigure4 figsize84 clearTrue pltsemilogyf SU pltgridTrue pltaxis0 12r 1e10 1e02 pltxlabelFrequência Hz pltylabelEspectro do deslocamento m²Hz sU2 nptrapz SU f sU4 nptrapzffSU f nu npsqrtsU4sU2 lnu npsqrt2nplog60nu Tempo de excitação é 60 segundos g lnu 05772lnu sU npsqrtsU2 up gsU print O fator de pico da resposta em deslocamento é 062f format g print Amplitude rms da resposta em deslocamento é 040f mm format 1000 print Valor de pico da resposta em deslocamento é 040f mm format 10 k 1000 uA Feq up k print Força estática equivalente é 041f N format Feq PEC00025 githubmmaiarochaPEC00025treemaster Classe16TestP2bipynb githubmmaiarochaPEC00025treemasterClass16TestP2bipynb Universidade Federal do Rio Grande do Sul UFRGS Programa de PósGraduação em Engenharia Civil PPGEC PEC00025 Introdução à Teoria da Vibração Teste P2 20211 múltiplos dof e sistemas contínuos NOME CARTÃO Instruções 1 Entregar a resolução da prova em arquivo único com no máximo 10Mb até às 12h de amanhã 01 de junho de 2021 2 Recomendase verificar atentamente se todas as folhas de resolução foram incluídas no arquivo gerado pois não serão aceitas entregas posteriores 3 Na primeira folha do arquivo deverá constar claramente o NOME e o cartão de MATRÍCULA 4 A consulta ao material de estudo e o uso do computador para cálculos são LIVRES 5 A prova deve ser realizada INDIVIDUALMENTE sem recorrer ao auxílio de colegas ou outras pessoas Caso se verifique o descumprimento desta regra todos os envolvidos na fraude terão uma nota da prova zerada Em 1 Importando módulos Python necessários para este notebook esta célula deve ser executada com shiftenter antes de qualquer outra célula P importar numpy como np importar scipylinalg como sc importar matplotlibpyplot como plt de MRPy importar Questão 1 Dados do problema Em 2 Função para cálculo dos modos de vibração Em 3 H 3 altura de cada pavimento m M 10000 massa de cada pavimento kg f1 1 frequência fundamental Hz zt 001 amortecimento modal adim mesma nos dois modos g 981 tração da gravidade ms2 def modosdevibração K M 1 Usa scipy para resolver o problema de autovalor padrão w2 Phi sc eig K M 2 Garanta a ordem crescente dos autovalores iw w2 argsort w2 w2 iw Phi Phi iw 3 Autovalores para frequências de vibração wk np sqrt np real w2 fk wk 2 np pi 4 Mass matrix normalization Mk npdiagnpdotPhiT npdotM Phi for k in rangelenwk Phik PhiknpsqrtMkk 5 Return results return fk wk Phi Monta matrizes e calcula modos In 4 Rigidez individual de cada barra 5168kNm Frequência no primeiro modo 100Hz Frequência no segundo modo 262Hz Visualiza modos In 5 K 100 rigidez de cada coluna incógnita KG Knparray2 2 2 4 rigidez global MG Mnparray1 0 0 1 massa global fk wk Phi vibrationmodesKG MG K f1fk02 determina a rigidez correta fk fknpsqrtK calcula todas as frequências wk fk2nppi em rads printRigidez individual de cada barra 061fkNmformatK1000 printFrequência no primeiro modo 062fHzformatfk0 printFrequência no segundo modo 062fHzformatfk1 printwk pltfigure1 figsize128 x Hnparange3 for k in range2 qk npzeros3 qk1 Phi1k qk npmaxnpabsqk adjust scale for unity amplitude pltsubplot12k1 pltplotqk x bo pltplotqk x pltxlim15 15 pltylabelstrk1 pltylim 00 70 plttitlefk 042fHzformatfkk pltgridTrue Questão 2 A excitação tem mesma amplitude nas frequências de 05 e 3Hz Para usar as amplificações dinâmicas vamos admitir que o pico das respostas modais poderão estar em fase In 6 Máximo deslocamento no pavimento superior 00398m O mesmo cálculo agora por simulação integrando por Fourier através do módulo MRPy In 7 Calcula forças modais e resolve equações de equilíbrio desacopladas In 8 Máximo deslocamento no pavimento superior 00388m FG 01gnpdiagMGreshape21 amplitude das forças nos Fk npmatmulPhiT FG amplitude das forças mod Mk npdiagnpdotPhiT npdotMG Phi massas modais Kk wkwkMk rigidezes modais uk npempty2 aloca memória para respo for k fn in enumeratefk bt 05 30fn frequências componentes AD npsqrt11 bt22 2ztbt2 respectivas amplificaçõe ukk FkkKkknpsumAD pico da resposta modal a u npmatmulPhiuk printMáximo deslocamento no pavimento superior 064fmformatu0 Td 32 N 1024 t nplinspace0 Td N domínio do tempo F FGnpsinnppit npsin6nppit força dinâmica Fk MRPynpmatmulPhiT F TdTd cria objeto MRPy for k in range2 Fkk Mkk prepara para solução uk FksdofFourierfk zt calcula respostas modais u npmatmulPhiuk deslocamento nos pavimen printMáximo deslocamento no pavimento superior 064fmformatu0max uplottime1 figsize105 A diferença dos dois resultados se deve a que o pico das respostas modais não está perfeitamente em fase Portanto a solução numérica que é a mais precisa apresenta amplitude ligeiramente menor Questão 3 Primeiro vamos calcular a resposta exata aplicando as condições de contorno na solução geral onde φx cos px cosh px cos px cosh px sin px sinh px si C1 C2 C3 C4 As condições de contorno são Aplicando essas condições na solução geral temos para Portanto e Por outro lado para Colocando as equações acima em forma matricial temos Fazendo o determinante da matrix de coeficientes igual a zero temos as frequência naturais Estas frequências podem ser calculadas numericamente como mostrado abaixo In 9 Cantilever beam frequency parameter is 1570796 Ou seja o parâmetro de frequência parece ser e a frequência fundamental resulta que coerentemente corresponde à frequência fundamental de uma viga biapoiada com vão Isso está correto já que a condição de apoio da direita equivale a uma condição de simetria para uma viga com o dobro do vão Agora vamos refazer o cálculo propondo a seguinte função aproximada para a forma modal p4 μ EI ω2 φ0 0 φ L φ L φ 0 0 0 0 x 0 φ0 0 φ 1 1 1 1 0 0 0 0 0 C1 C2 C3 C4 1 1 1 1 0 0 0 0 0 C1 C2 C3 C4 C1 0 C2 0 x L L φ L φ cos pL cosh pL cos pL cosh pL 0 C3 C4 cos pL cosh pL cos pL cosh pL 0 C3 C4 cos pL cosh pL cos pL cosh pL cos pL cosh pL cos pL cosh pL C3 C4 0 0 π2 ω1 π 2L 2 EI μ 2L φx 2Lx 1 L2 x2 def chareqx x x0 A nparray npcosxnpcoshx npcosxnpcoshx npcosxnpcoshx npcosxnpcoshx return nplinalgdetA from scipyoptimize import fsolve p fsolvechareq 10 printCantilever beam frequency parameter is 086fformatp0 ou seja uma parábola que apresenta derivada nula para e portanto respeita algumas condições de contorno A escala desta forma modal é intencionalmente escolhida como sendo unitária As derivadas dessa forma modal aproximada são Observase que a função proposta não cumpre a condição de momento nulo na extremidade da esquerda mas vamos em frente A correspondente energia cinética de referência é enquanto a energia potencial elástica resulta Portanto o quociente de Rayleigh resulta onde e portanto esse valor apresenta um erro de aproximadamente 54 esperado em relação ao valor correto e um erro de aproximadamente 11 em relação à frequência correta Observe que o quociente de Rayleigh sempre fornece frequência igual ou maior que a exata Questão 4 Dados do problema x L x φ x φ 2L 2x L2 2 L2 μ x dx μL Tref 1 2 L 0 φ2 4 15 V EI dx 1 2 L 0 φx 2 2EI L3 ω1 V Tref 2EI 15 4μL4 12014 2L 2 EI μ 331 12014 π In 10 Vamos considerar a resposta apenas no primeiro modo A dissipação de energia por amortecimento é desprezada e a energia total do sistema deve se manter constante e igual à energia potencial gravitacional da pessoa no início da queda Por questão de simplicidade admita que a viga já está deformada pelo peso próprio quando se determina a altura de queda da pessoa Também vamos considerar que o choque é perfeitamente inelástico ou seja a viga e a pessoa unida seguem após o contato Observe que a forma proposta modal é normalizada pela unidade de modo que ela tem valor unitário na extremidade da direita Desta forma deslocamento vertical e deslocamento modal tem mesmo valor numérico no ponto B A energia cinética total do sistema após o choque é calculada como onde é a velocidade inicial da extremidade direita da viga logo após o choque que é numericamente igual à velocidade inicial no espaço modal Igualandose as energias chegase a A frequência natural no primeiro modo precisa ser recalculada pois agora a viga tem também a massa da pessoa interna na extremidade direita A nova energia cinética de referência é A energia potencial hidrogênio permanece a mesma e portanto a frequência natural resultado menor Sem a massa da pessoa interna a frequência natural calculada na seção anterior seria de A amplitude total do deslocamento modal é a soma da amplitude devida à velocidade inicial com o deslocamento devido à carga impulsiva Dada a escala unitária da forma modal a força modal tem o mesmo módulo de força aplicada na extremidade da direita O formato retangular da carga impulsiva choque inelástico implica que o fator de amplificação dinâmica da resposta estática é igual a 2 Para calcular a resposta estática é necessário conhecer a massa modal E mgh 7848J φL 1 T μ dx m 360 1 2 L 0 φx v0 2 1 2 φL v0 2 v2 0 v0 E T 148ms v0 7848 360 μ x dx m L 360 Tref 1 2 L 0 φ2 1 2 φ2 V 3042rads ωn V Tref 2 36 106 360 63 3227rads A uBest A uBmax v0 ωn uBest L 6 comprimento da viga m m 80 massa da pessoa kg mu 200 massa por unidade de comprimento kgm EI 36e6 estreitamento à flexão Nm2 Lembrando que a dor modal é dada por a resposta estática é calculada como Substituindo valores Em M μ x dx m L 720kg L 0 φ2 φ2 K M ω2n 118mm uBest mgφL K 80 981 1 720 30422 2 000118 51cm uBmax 148 3042 PEC00025 githubmmaiarochaPEC00025treemaster Classe16TestP2aipynb githubmmaiarochaPEC00025treemasterClass16TestP2aipynb Universidade Federal do Rio Grande do Sul UFRGS Programa de PósGraduação em Engenharia Civil PPGEC PEC00025 Introdução à Teoria da Vibração Aula 16 Teste P2 múltiplos graus de liberdade e sistemas contínuos P22019 Pergunta 1 Pergunta 2 Pergunta 3 Pergunta 4 Prof Marcelo M Rocha Drtechn ORCID httpsorcidorg0000000156401020 Porto Alegre RS Brasil Em 1 P22019 Observação esta prova deve ser resolvida com o auxílio de uma calculadora científica que deve ser capaz de resolver autoproblemas sistemas lineares e integrais O tempo total disponível para a resolução da prova é de 2h duas horas O aluno poderá preparar uma folha A4 frente e verso com as informações a serem consultadas durante a prova Questão 1 Um sistema estrutural é modelado como um sistema discreto de dois graus de liberdade conforme mostrado na figura Cada coluna possui rigidez à flexãoEEU 500kNm 2 comprimentoeu 4m e são assumidas como não tendo massa relevante As vigas do piso são assumidas como perfeitamente rígidas e com massa total aglomeradam 4toncada Supõese que o sistema apresente um amortecimento viscoso com razão deζ 001em todos os modos de vibração 1 Defina as matrizes de rigidez massa e sistema de amortecimento 1 pts 2 Determine e esboce os dois modos de vibração naturais indicando as frequências de vibração associadas 2 pts Importando módulos Python necessários para este notebook esta célula deve ser executada com shiftenter antes de qualquer outra célula P importar numpy como np importar matplotlibpyplot como plt importar scipylinalg como sc Resposta As matrizes de rigidez e massa são Em 2 Global stiffness matrix 187500 187500 187500 375000 Global mass matrix 4000 0 0 4000 To specify the damping matrix we must first calculate the vibration modes and frequencies EI 500000 rigidez à flexão de uma única coluna Nm2 m 4000 massa de um único piso kg L 4 comprimento da coluna m k 12 EI L 3 rigidez de coluna única KG np array 2 k 2 k 2 k 4 k matriz de rigidez global MG np array m 0 0 m matriz de massa global printGlobal stiffness matrix KG print Global mass matrix MG In 3 And now we can calculate the coefficients that multiply the stiffness and mass matrices to build a Rayleigh damping matrix that is also orthogonalized by the eigenvectors Uses scipy to solve the standard eigenvalue problem w2 Phi sceigKG MG Ensure ascending order of eigenvalues iw w2argsort w2 w2iw Phi Phiiw Eigenvalues to vibration frequencies wk npsqrtnprealw2 fk wk2nppi pltfigure1 figsize86 clearTrue x nparange0124 for k in range2 pk npzeros3 pk1 Phi1k pk npmaxnpabspk adjust scale for unity amplitude pltsubplot12k1 pltplotpk x bo pltplotpk x pltxlim15 15 pltylabelstrk1 pltylim 00 10 plttitlefk 042fHzformatfkk pltgridTrue In 4 Mass matrix coefficient a0 006124 Stiffness matrix coefficient a1 000131 Rayleigh damping matrix original 48989794856 24494897428 24494897428 73484692283 Rayleigh damping matrix orthogonalized 33851115694 0 0 88623371445 Question 2 The system is now subjected to an initial kinematic condition which consists of an imposed displacement on the lower floor u20 1cm only and it is then released to vibrate Accounting for the two vibration modes calculate the peak displacement and the peak acceleration at the system upper floor caused by this initial condition 2 pts Answer For the modal superposition we must firstly calculate the modal masses and the modal stiffnesses In 5 Modal masses 4000 4000 Modal stiffnesses 71619 490881 The initial condition is of displacement type no initial velocity what implies a cosine type response Recalling that Φ is a orthogonal matrix it means that its transpose is equal to its inverse ut Φ ukt ukt Φ ut where ut is the nodal response and ukt is the modal response The initial modal displacements are simply given by zeta nparray001 001 damping for two modes i and j A nparray1wk0 wk0 1wk1 wk12 alpha nplinalgsolveA zeta CG alpha0MG alpha1KG Rayleigh viscous damping matrix printMass matrix coefficient a0 065fformatalpha0 printStiffness matrix coefficient a1 065fformatalpha1 print Rayleigh damping matrix original CG print Rayleigh damping matrix orthogonalized npdotPhiT npdotCG Phi Kk npdiagnpdotPhiT npdotKG Phi Mk npdiagnpdotPhiT npdotMG Phi printModal masses 060f 160fformatMk printModal stiffnesses 060f 160fformatKk In 6 Initial modal displacement at mode 1 0005257 Initial modal displacement at mode 2 0008507 Initial modal displacement at mode 1 0005257 Initial modal displacement at mode 2 0008507 The total response is a superposition of modal responses which are cosine functions with the respective frequencies and amplitudes In 7 The accelerations are obtained from the twofold derivative of the cosine sum u0 nparray000 001 initial displacements in nodal coordinates u0k npdotPhiT u0 initial displacements in modal coordinates printInitial modal displacement at mode 1 086fformatu0k0 printInitial modal displacement at mode 2 086fformatu0k1 The most general way as shown in classroom considering phase is pi2 u01 npdotPhi0 npdotMG u0npsinnppi2Mk0 u02 npdotPhi1 npdotMG u0npsinnppi2Mk1 printInitial modal displacement at mode 1 086fformatu01 printInitial modal displacement at mode 2 086fformatu02 T 10 N 200 t nplinspace0 T N time domain uk nparrayu0k0npcoswk0t u0k1npcoswk1t modal responses u npdotPhi uk100 total responses cm pltfigure2 figsize12 4 clearTrue pltplott u0 b t u1 r pltxlim 0 T pltxlabeltime s pltylim2 2 pltylabelut cm pltlegendupperlower pltgridTrue In 8 Finnaly answering the question the peak displacement and acceleration amplitudes in the upper floor are In 9 Peak upper displacement 0876cm Peak upper acceleration 0064G It can be seen that as expected the second mode dominate the structural response Question 3 The cantilever beam shown in the figure has a constant flexural stiffness EI 1000kNm2 and mass per unit length μ 200kgm 1 Propose a function that resembles the first vibration mode Calculate the associated potential elastic energy V and the reference kinetic energy Tref With these energies estimate the natural vibration frequency for the first mode using the Rayleigh quocient 2 pts 2 Calculate the modal mass and the modal stiffness and then use these parameters to estimate the static displacement at the cantilever tips caused by a point load W 10kN placed at this same position 1 pts ak nparrayu0k0wk0wk0npcoswk0t u0k1wk1wk1npcoswk1t modal accelerations a npdotPhi ak981 nodal accelerations G pltfigure3 figsize12 4 clearTrue pltplott a0 b t a1 r pltxlim 0 T pltxlabeltime s pltylim02 02 pltylabelat G pltlegendupperlower pltgridTrue printPeak upper displacement 053fcmformatu0max printPeak upper acceleration 053fG formata0max Answer We will try and compare two different solutions a parabolic and a sinusoidal functions They are φ1x 1 27 x 3x 9 and φ2x 1 2 sin πx 12 Both tentative solutions respect the kinetic condition of zero displacement at supports located ate coordinates x 3m and x 9m The script below shows a comparison plot In 10 EI 1000000 flexural stiffness mu 200 mass per unit length L 12 total length N 200 number of segments X nplinspace0 L N length discretization ph1 lambda x x 3x 927 first solution ph2 lambda x 1 npsqrt2npsinnppix12 second solution pltfigure4 figsize12 4 clearTrue pltplotX ph1X b X ph2X r pltplot3 9 0 0 bo pltxlim 0 L pltxlabelx m pltylim2 2 pltylabelphi nondim pltlegendphi1phi2 pltgridTrue The sine function has an important feature which is zero curvature at cantilever tips where bending moments must be zero The parabolic function is the simplest but presents constant curvature along all beam length The rotations are calculated as ϕ 1x 1 272x 12 and ϕ 2x π2 12 cos πx 12 while the curvatures are given by ϕ 1x 2 27 and ϕ 2x π22 144 sin πx 12 The script below compares the curvatures for each solution In 11 ph1xx lambda x 227x0 first solution ph2xx lambda x nppinppinpsqrt2144npsinnppix12 second solutio pltfigure5 figsize12 4 clearTrue pltplotX ph1xxX b X ph2xxX r pltxlim 0 L pltxlabelx m pltylim005 015 pltylabelphixx 1m2 pltlegendphi1phi2 pltgridTrue The curvatures are quite close at the center but it is overestimated by the parabolic function at the cantilever tips The potential elastic and the reference kinetic energy finally are In 12 Potential elastic energy for solution 1 327572J Potential elastic energy for solution 2 280446J Reference kinetic energy for solution 1 2035J Reference kinetic energy for solution 2 2381J And the natural vibration frequencies estimated with Rayleigh quotient are In 13 Natural vibration frequency for solution 1 202Hz Natural vibration frequency for solution 2 173Hz If one recalls that the true vibration mode minimizes the Rayleigh quotient the lowest value obtained with the sinusoidal function is likely to be closer to the exact solution The relative error between both tentative functions is approximately 17 and the correct natural frequency must be a little below 173Hz Now we will proceed with the calculation of modal mass and modal stiffness dx LN V1 EInptrapzph1xxXph1xxX dxdx2 V2 EInptrapzph2xxXph2xxX dxdx2 T1 munptrapz ph1Xph1X dxdx2 T2 munptrapz ph2Xph2X dxdx2 printPotential elastic energy for solution 1 051fJformatV1 printPotential elastic energy for solution 2 051fJ formatV2 printReference kinetic energy for solution 1 051fJformatT1 printReference kinetic energy for solution 2 051fJformatT2 wn1 npsqrtV1T1 wn2 npsqrtV2T2 fn1 wn12nppi fn2 wn22nppi printNatural vibration frequency for solution 1 052fHzformatfn1 printNatural vibration frequency for solution 2 052fHzformatfn2 In 14 For static analysis the modal displacement is obtained from modal force divided by modal stiffness In 15 Static displacement of cantilever tip for solution 1 509cm Static displacement of cantilever tip for solution 2 738cm The error in the displacement at cantilever tip for the two solutions is quite high over 40 for the two tentative functions diverge noticeably in that position we recommend this result to checked with the Ftool software A comparison of displacement solutions for the whole beam is shown below In 16 Mk1 munptrapzph1Xph1X dxdx modal mass and Kk1 Mk1wn12 stiffness for solution 1 Mk2 munptrapzph2Xph2X dxdx modal mass and Kk2 Mk2wn22 stiffness for solution 2 W 10000 point load downwards Fk1 Wph16 modal static force Fk2 Wph26 uk1 Fk1Kk1 modal displacement uk2 Fk2Kk2 u1 uk1ph112 displacement at cantilever tip u2 uk2ph212 printStatic displacement of cantilever tip for solution 1 052fcmformatu11 printStatic displacement of cantilever tip for solution 2 052fcmformatu21 pltfigure6 figsize12 4 clearTrue pltplotX 100uk1ph1X b X 100uk2ph2X r pltplot3 9 0 0 bo pltxlim 0 L pltxlabelx m pltylim6 12 pltylabelu cm pltlegendphi1phi2 pltgridTrue Question 4 The same point load from previous question is now applied suddenly from zero to its final magnitude what causes a dynamic amplification on the beam displacements Estimate the peak displacement and the peak acceleration at the cantilever tip 2 pts Answer The solution for some impulsive loading is well known to be the static solution multiplied by a dynamic amplification factor In the case of a step load Heavisides function this amplification factor is 2 Hence In 17 Deslocamento dinâmico da ponta do cantilever para solução 1 1018 cm Deslocamento dinâmico da ponta do cantilever para solução 2 1477 cm As acelerações de pico são Em 18 Aceleração na ponta do cantilever para solução 1 0835G Aceleração na ponta do cantilever para a solução 2 0887G Observase que o erro na resposta de aceleração na ponta do cantilever de aproximadamente 6 não é tão alto quanto o da resposta de deslocamento print Deslocamento dinâmico da ponta do cantilever para solução 1 052f cm print Deslocamento dinâmico da ponta do cantilever para solução 2 052f cm ak1 uk1 wn1 2 ak2 uk2 wn2 2 a1 ak1 ph1 12 a2 ak2 ph2 12 print Aceleração na ponta do cantilever para solução 1 053f G format a1 print Aceleração na ponta do cantilever para solução 2 053f G format a2