Transcription
AnáliseNuméricado ModeloparaVibraçõesde MembranasFronteiraRaquel Defelippo Rodriguesde Kirchhoff-CarrierElásticascomMóveleMauro Antonio RinconResumoNeste trabalho, estudamos um modelo de vibrações transversais emmembranas elásticas que foi desenvolvido por Medeiros e FerreI como umaextensão do modelo de Kirchhoff-Carrier, uma vez que leva em contaa troca de tamanho da membrana durante a vibração. O modelo deMedeiros e Ferre] não possui solução analítica, então usamos o Métododos Elementos Finitos para a resolução do modelo no espaçoe um Métodode Diferenças Central, o Método de Newmark, para a resolução do modelono tempo e assim obtemos a solução aproximada do problema. Algunsexemplos numéricos são apresentados.Palavras Chave: Equação de Kirchhoff-Carrier, extremidades móveis,solução numérica1IntroduçãoA descrição matemática para pequenas vibrações transversais de cordas elásticascom extremidades fixas é uma velha questão.Os primeiros resultados matemáticos para o estudo de pequenas vibraçõestransversais em cordas elásticas foram obtidos por Taylor em 1713 e informaçõesmais precisas datam de 1743 por D IAlembert e posteriormente por Euler .O modelo proposto por D' Alembert é: C2 ât2ondec2 !.- I sendomT a tensão(1)âX2constantedurantea vibraçãoe m a massadacorda.Para obtermos o modelo de D' Alembert é preciso impor algumas restriçõesao problema físico.1
Outros modelos de problemas físicos de pequenos deslocamentos verticaisem cordas elásticas com extremidades fixas foram propostos por Kirchhoff eCarrier. Sendo o modelo de Kirchhoff-Carrier dado por:82u""ã"i2- (LTOk [8u(ã;))2dx82uã;2 O(2)É importante observar que o modelo de D' Alembert não considera os pequenosdeslocamentoshorizontais aqui considerado. Essepequeno deslocamentohorizontal, no modelo matemático, é representadopelo termo não linear .Em 1999, Medeiros e FerreI propuseram em [3] um modelo que é uma extensão do modelo de Kirchhoff-Carrier, pois considera a troca de tamanho dascordas elásticas durante a vibração apresentadopor: - [ '"Y(t) -Lo {3(t)!8t2mmLo( )8x2m'Y(t)]2dx 8X2 O(3)a(t)SendoLo o tamanho inicial, TOa tensão inicial, m a massa,k o módulo de Youngdo material vezesa seçãotransversal da corda e u(x, t) o deslocamentoverticaldo ponto x da corda no tempo t. As funções a(t) e [3(t) dão o movimento dosextremos da corda durante a vibração sendo a(t) [3(t) e '"Y(t) [3(t) -a(t).Em 2000, Medeiros e FerreI em [4] extenderam o modelo (3) para o casobidimensional das membranas elásticas com sua fronteira móvel. O modelopara o caso das membranas elásticas é dado por:[ TO k K2(t) -K 82U(:J;,8t2 t) - K kr2] 2m1TK2(t) Jn. IVu(:J;,t) I dx 6.u(x, t) 0(4)Sendo nt um domínio circular de raio K(t), TOa tensão inicial, m a massa,Ek 2(i"" ãj" onde E é o módulo de Young do material e u o coeficiente dePoisson. O deslocamentovertical do ponto :J; (Xl,X2) E 1R2da membrana notempo t é representadopor u(:J;,t).2Formulaçãodoproblemao modelo para pequenas vibrações transversais de membranas elásticas comfronteira livre foi apresentadocomo -[a(t) b(t) k.lvuI2dX]6.u O(5)onde as funções a(t) e b(t) são dadas pora(t) mmK2o2(K2(t)-K )(6)
kb(t) 2m7rK2(t)(7)sendo que a(t), b(t) e o domínio de integração nt , um disco com raio K(t),variáveis com o tempo t. Onde u u(:1:, t) representa o pequeno deslocamentovertical do ponto :I: (Xl,X2) da membrana nt.Consideremos uma membrana elástica n E JR2 de raio unitário e centro naorigem, isto én {(Xl,X2) E JR2;x x 1}(8)o domínio de integração nt representa a deformação do disco unitáriodada pela função K(t), tal que nt é definido 'v't O como sendo:nt {(Xl,X2)E JR2;:I: K(t)y,'v'y (Yl,Y2) E n}n(9)Observe que nt, para cada t ?: O, é um disco de raio K(t) e centro na origem.Representamos por no o disco Kon onde Ko K(O).Consideremos o domínio não cilíndrico Q do JR2 definido porQ Unt x {t}(10)rt x {t}(11)O t Te sua fronteira lateral Ê definida porÊ UO t Tonde rt denota a fronteira de nt. Supondo que rt e Ê são regulares.Tomando-se o operador diferencial definido pela função real u(:1:, t)iu(:1:,t) -[a(t) b(t) h,I'\lu (:I:, t)12dX] Âu(:1:, t)(12)onde a(t) e b(t) são dados por (6) e (7), respectivamente.O problema que estudaremos será o de determinar uma função u(:1:,t), noespaço das soluções HJ(nJ n H2(nt), tal que{I iu(:1:,t)u(:1:,t) f(:1:,t),u(:1:,O) uo(:1:),as equações (13) formam(:1:,t) E Q 0,(:1:,t) E Êâu8t"(:1:,O)o que chamaremos3 Ul(:1:),no Konde Problema(I), onde:l:(13) (Xl, X2).
3ProblemaequivalenteTrabalhar diretamente com o operador Lu(x, t), apresentado em (12), originacertas dificuldades devido a natureza variável do domínio de integração nt,dado em (9), com relação ao tempo t, o que torna o Problema (I), (13), difícilde ser solucionado diretamente.Passaremos então a transformar o domínionão cilíndrico Q num domínio cilíndrico Q n x (O,T), onde S1 é o discounitário de centro na origem definido em (8). Isso é feito por meio da seguintetransformação:A:ondeQ(x, t)-tH-Q(y, t)(14)xy K(t)(15)e K : [0, T) -t 1R é uma função de classe C2 ([0, T); 1R) satisfazendo as seguintescondições.Assume-se que K(t)satisfaz as condições:K E C2 ([O,T[j IR) ,K(t)I K"(t)e O K'(t)?: Ko ?: 1I Co Kl para todo 0 t T. Onde Kl satisfaz:K2 12medT(CO,C1,C2,C3,C4,Ko,K1).com 4( TO -k )TO k.Observação:As constantesC1, C2, C3, C4 são definidasKO 211vol12m OCl IVll2 como 211vO114jm1l"K OIllpIIH2(O) C2IdlpIL2(O)IIvIIHJ(O) C3IdvIL2(O)[1 12dr C41dulO parâmetro d 0. A função T é definida por:d2 (ao 1)204 2T 4 [C3(C2 ao 2) 2 ao] Kl 4KoKl 2C (1 ao) K34(C3 C2)2 (1 CO)2(C3 02)2K21 dK2OO -42m4
Maiores detalhes em [4].Aplicando a transformação A definida em (14) no operador Lu(x, t), dadoem (12), obtemos o operador Lv(y, t).Lv(y,t) a2v-a2ta2va2va a-ã(t)6.v bij(y,t)-Yi Ci(y,t)-Yjavata ei(Y,t)- aYiYi(16)ondeã(t) bij(Y, t) ( K(i) ) yiYjCi(Y, t) -2ei(Y, t) rl2A transformação(a(t) b(t)iIV'VrdY)(17)(18)KI(t»( K(i)KI(t» ) Yi( KI(t» ) 2K(i)(19)K1'(t) 1-K(t)J yi(20)A dada pela seguinte proposição.Proposição1: O operador Lv(y, t), definido em (16), é obtido a partir doopemdor Lu(x, t), definido em (12), por meio da mudança de variável dada pelatransformação(14).Obtemos assim, em um domínio retangular Q n x (0, r),equivalente de determinar no espaço das soluções HJ(n) nH2(n),.v v(y,t), com y (YI,y2), tal que,{II Lv(y,t) g(y,t),v(y,t) O,av(y,O) vo(y),-ai(y,o)(y,t)EQyEÊ VI (y),o problemauma função(21)y Enas equações (21) formam o que chamaremos de Problema equivalente (II).Em [4], provou-se a existência e unicidade de solução do Problema equivalente (II). Como a transformação A, dada em (14), é um isomorfismo então asolução v(y, t) do Problema (II) é transformda na solução u(x, t) do Problema(I).Por conveniência a análise numérica usando elementos finitos no espaçoe diferenças finitas no tempo será baseada no Problema equivalente (II) dedomínio cilíndrico, ou seja, determinaremos a solução aproximada do Problemaequivalente (II).5
4MétododosElementosFinitosFormulação VariacionalComo o Método dos elementos finitos não é aplicável diretamente ao Problema (II), (21), passaremosa construção da formulação variacional.Seja D(n) o espaço das funções teste, infinitamente diferenciáveis com suporte compácto em n e w E D(n). Multiplicando a primeira equação do Problema (II), (21), por w e integrando emn {(Yl,Y2)/Y Y 1}obtemos:{ 82V{{82vJn ai2wdy -Jn ã(t)Âvwdy Jn bij(Y,t) wdy{82V{8v{ Jn Ci(Y,t)8"ta"Y:wdy Jn ei(Y, t)ay;wdy Jn g(y, t)wdy(22)Método de GalerkinConsideremos Vm C HJ(n)nH2(n)um subespaçode dimensão finita gerado pelos m primeiros elementosda base do espaçoHJ(n) n H2(n), ou sejaVm [ Pl, P2, P3,.", pm]onde [ Pi]iEJVé uma base de HJ(n) n H2(n), Agora passamosa buscar umasolução aproximada vh(y, t) do Problema (II) no subespaçoVm.Problema AproximadoAproximamos o Problema (II), (21), pelo problema de determinar no espaçodas soluçõesVm, uma função v h vh(y,t), com y (Yl,Y2), tal que,{III Lvh(y, t) g(y, t)hv(y, t) E QA(y,t) Ov(y,O) vo(y),yEE8va-i(Y'O) v 1(y),(23)YEnAs equações (23) formam o que chamaremos de Problema (III).Substituindo v h vh(y, t) E Vm em (16) temos o operador Lvh(y, t) dadopor:Lvh(y,t) 82vh-8t282vh-ã(t)ÂVh bij(y,t)-82Vh8Yi8Yj Ci(Y't)- 8 t 8Yi ei(Y,t)-8vhaYi (24)onde vh(y, t) E Vm é uma soluçãoaproximada de v(y, t), e ã(t), bij (y, t), Ci(Y,t)e ei(Y, t) são dados por (17), (18), (19) e (20), respectivamente.6
Assim, (22) é agora dada porla-a2vh2nlwdy-tã(t) vhwdy1 n 1 Ci(Y,t)-a2vhat ayinbij(Y,t)-a2vhayi aYjn1 ei(Y,t)-wdy avhawdyYin1 wdyg(y,t)wdy(25)nAplicando o teorema da divergência na segunda e na terceira integrais,lembrando que bij(y,t) é dado em (18), obtemosJn wdy ã(t) Jn VvhVwdyr a2vh( KI(t) ) 2-K(i)ravh( K(i)KI(t)) 2 Jnr Yi-av;-wdyavh-3rawera2vhJnYiYj-av;-8Yjdy JnCi(Y,t) wdy l ei(Y,t)-avhanwdy1 yig(y,t)wdy(26)nCom vh(y, t) E Vm podemos escrever vh(y, t) da seguinte maneiramvh(y, t) L dk(t) pk (y),k l Pk(Y) E Vm(27)Para obtermos a solução aproximada vh(y, t) E Vm é necessário determinaros coeficientes dk(t).Substituindo (27) em (26) obtemosi-t( dk Pk(Y))WdY ã(t) i)3( 21(t )(K K(t )) ka.dw.Jn ny .1.ak l) wdy (dk k lyY.) 1( ) (d Y.Y, L., k a r Ci(Y,t) fY.) dy .yr ei(Y, t) fJnk l(dk Y.)g(y, t)wdyObservemosobtendo-seque a funçãog(y,t)podeser interpoladam Lgk(t) Pk(Y)k lt) gk(t),wdy(28)g(y,t)onde g(yk'( dk V pk(Y) ) Vwdy2- i(d. Y. L.,nk ltpara k 1, 2, ., m7pela função po(Y)
Além do mais, tomando em particular w P1(Y)E Vm em (28) obtemos{ dk (i Pk(Y) PI(Y)dy) ã(t)dk (i( KI(t) )2-3K"(i)() (dk2 dkrJn(rÔ pk(Y))Yi PI(Y)dy(YiYj) dy-K(t)JnÔYi dk (iCi(Y,t) PI(Y)dY) i( 9k(t) pk(Y) ) P1(Y)dy V Pk(Y)VCPI(Y)dy)ôYj) dk (iei(Y,t) cPl(Y)dY)}(29)para l 1,2, .m5DiscretizaçãodoDomínioConsideremos uma partição do domínio !1 em subregiões !1e, de tal forma asatisfazer as condições:Nel!1 U !1e,e l!1e n !1k 0, e konde N el é o número total de elementos. Na partição do domínio definimos.os nós globais A, A 1, 2, ., Nnos, onde Nnos é o número total de nós damalha. Tal partição é necessária para que possamos encontrar a solução vh(y, t)do Problema (III).A malha utilizada está apresentada na Figura 1, onde está apresentada umamembrana circular dividida em 73 nós e 128 elementosOptamos por trabalhar com subregiões !1e triangulares devido a geometriado domínio, considerando assim 3 nós locais em cada elemento finito ne.Como estabelecemos um mapeamento biunívoco entre a numeração globale a numeração local dos nós é possível obter a solução em todo n a partir dassoluções locais em !1e.Para cada elemento finito !1e(y), define-se as funções de interpolação local P: com as seguintes propriedades:{ p:(y) 1,y y: E neO,0,y ygE!1e,y lt !1e8b a(30)
o.00.I.I00.5 0.5Figure 1: Membrana dividida em 73 nós e 128 elementosonde y e yb são as coordenadas locais de ne. Pela aplicação do Método deGalerkin, vimos em {27) que a solução aproximada vh{y, t) é dada pormVh{y, t) L dk{t)'Pk{Y),k l'v'y E nLocalmente a solução aproximada é dada porv {y,t)3 L d {t)'P {Y),a lDesta forma a solução aproximada Vh{y,t)por'v'Y E nedo Problema {III),{31){23), é dadaNel U v (y,t)(32)e londe N el m é o número de elementos finitos.Considerando a notação acima passamos a definir {29) em termos de cadaelementovh{y,t)9
{ d: (ke CP:(Y)CPb(Y)dY) ã(t) d: (ke VcP:(Y)VcPb(Y)dy)-3 ( )2d:-( K(t) ) d: (k.(k. Yi CPb(Y)dY)(Jn.r (YI.YJ.) d8Yi2deoY)8YjCi(y,t) cpb(Y)dY) d: (k.ei(Y,t) lpb(Y)dY)}3 Li0 1g:(t)Ip:(Y)Cpb(Y)dy(33)n.para b 1, 2, 3 e e 1, 2, ., N el.A equação (33) é de difícil solução devido a existência das funções Ci(Y,t),ei(Y, t), yi e Yj nas integrais.Com o conceito local apresentado anteriormente e assumindo que a área decada elemento ne é suficientemente pequena passaremos a considerar as funçõesCi(Y, t), ei(Y, t), Yi e Yj constantes em cada elemento finito ne.Para isso, consideramosYie 31 (Yi1 Yi2 Yi3 ) e( 34)para 1 i, j 2Assim,temos { d: k. cP:(y)cpb(y)dy d: ã(t) k. VcP:(y)\7cpb(Y)dy-3 ( ) 2yf k. cpb(Y)dy() -deo d:2( ': I!) rYI YJ Jn.K(t)ci(t)3 Lg:(t)0 1rJn.1 d8Yiy(.Y)cpb(Y)dY d:I8Yjyei(t) r Jn.cP:(y)cpb(y)dyn.para b 1,2,3 e e 1,2,.,Nel.Colocando em evidência os termos comuns, obtemos:10y(.Y)lpb(Y)dyI}
1 P (Y) Pb(Y)dy 3{ d [LO. d [ã(t)L.V p (y)V pb(Y)dy-3( K(t) ) 2(y y ).1-[ c (t)]a l3 { g:(t) U.rJo.r pb(Y)dyJo.Y.( )2YiL. dy e (t)âyi]âyj.Jo. pb(Y)dYr .(y)dâyi Pby]} P (Y)'Pb(Y)dY)}(35)para b 1,2,3 e e 1,2,.,Nel.O queé equivalenteat{ d [L. P:(Y) Pb(Y)dy] d [-2 ( )YiL. 'Pb(Y)dy] d [ã(t) L. V p (y)V pb(Y)dy-3 ( )2-( K(t) ) 2(yiY )1((2) t)2 - K(t) r dyJo.âyi Pb(Y)dyây;y rK(t)yiL. P.(Y)dy.Jo.âyib]}{ g:(t) U. 'P (Y) Pb(Y)dY)}parab 1,2,3e e 1,2,.,Nel.Logo,t{d: [L. p:(y) pb(Y)dy] [-2( )yiL.( [ ã(t) ( -( rJo.V p:(y)V pb(Y)dy) 2 - 3 { g:(t) U.-) 2 (yiy ) r dyK(t)) yi L. Pb(Y)dY] P:(Y) Pb(Y)dY)}parab 1,2,3ee 1,2,.,Nel.11 Pb(Y)dy]1 Jo.âyiây;}(36)
Denotando por M, a matriz de ordem ( N el x N el) denominada matriz massa,a, a matriz de ordem (Nel x Nel) denominada matriz de amortecimento e K,a matriz de ordem (N el x N el) denominada matriz de rigidez e a o vetor deordem ( N el x 1) denominado vetor força. DefinimosNelM L M:b,e lNela L a:b,e lNelK L K:be lNele a L a e londe M:b, a:b e K:b são matrizes de ordem (3 x 3), e o vetor a:(3 x 1).As matrizes M:b, C:b e K:b extraídas de (36) são dadas por:M:b (37)é de ordemr (y) b(y)dyJn.a:b c (t) r (38)y (.Y) b(y)dyIJn.(39)K:b ( ã(t) h. V (y)V g(y)dy)-(\( K(t) ) 2(y y )r dy3 Jn. ayi (( -\\K:b \)ayj( ) 2- \ y r e(y)dyK(t)K(t) ) I Jn. ayi b(ã(t)r)V (y)V b(Y)dY-(b j(t)Jn.r\) dyJn.YIY3 ( he(t) h. b(Y)dY)(40)onde ã(t) é dada em (17)Denotandobe(t)-2b j(t). 1,3 1K(t)KI(t)() ( ) 2 ( Y Y2)2(t)( --e- ji ( K(t)KI(t)) 2 . 2 -2je(t) r (y yj)1,3 12(YIYl 2 YIY2 Y2Y2)e(l-eKI(t»e)K(t)) 2 - lK(t)K(t)12e ee(41)2K"(t)1-KWJrl)J(ye ye12e Yi2)(42)
Assim, K:b é dada porK:b (ã(t)fV p:(Y)V Pb(Y)dY) - (be(t)Jn. (fe(t)fJn.f dY8Yi8Yj) Pg(y)dY)10.Restringindo(43)8Yio vetor força a a cada elemento finito {1e, temosa: g:(t)f cp:(y) pb(Y)dyJn.(44)para 1 a, b 36FunçãodeInterpolaçãoConsideremos um triângulo {1e C {1 representado na Figura 2.eeYJY;Figure 2: Elemento triangular{1ePara cada {1e definiremos as funções de interpolação cP: contínuas no interiorde cada elemento de fronteira re, onde a 1,2,3 são os nós de cada triânguloe e 1, 2, ., N el são os elementos da malha.Para cada lado do triângulo definimos uma função de interpolação P: ondeas seguintes condições serão impostas:e)e(Yb Po {1,O,a ba b(45 )onde 1 a, b 3 e yg «yl)b, (y2)b).Observação:Para simplificar , passaremos a utilizar a seguinte notaçãoY (x,y) {:} Y (Yl,Y2)13(46)
É importante fazer a distinção entre as soluções u(:1:,t), onde :I: (Xl,X2)v(y,t), onde y (Yl,Y2) (x,y) dos Problemas (I) e (II), respectivamente.eAssim, de (46) em (45), temos que yg (xg,yg).Queremos determinar tp (y), tp (y) e tp (y), interpoladores satisfazendo (45).Em particular , quando fie é O triângulo retângulo apresentado na Figura 3, afunção de interpolação tPa(X,y) pode ser definida como a seguinte função lineartp (x,y) 1-x-y(47)tp2(x,y) x(48)tp3(x,y) y(49)onde tPl (y) , tP2(y) e tP3(y) são obtidas nas respecticas coordenadas:Yl (0,0),Y2 (1,0),Y3 (0, 1)y; .y; y;Figure 3: Elemento triangularfieComo tp (x,y) é uma função linear por partes esta pode ser representadacomo na Figura 4.14
rP'y;Figure 4: Representaçãogeométrica cP:(X, y)Podemos observar que o polinômio interpolador CPa(x,y),a 1,2,3 no elemento padrão fib é bem mais simples, assim é convenientefazer uma parametrização entre os diferentes elementosfie para um único elemento fib, onde por fibestamos denotando o triângulo retângulo apresentadona Figura 5.(0,1) (0.0)(1,0)Figure 5: Elemento triângulo padrão fib7TransformaçãoIsoparamétricaTrabalhar diretamente com os elementos triangulares fie tornaria o processomuito complicado, uma vez que as funções de interpolação deveriam ser calculadas para cada elemento. Assim, é conveniente trabalhar com elementostriangulares padrões fib, que foram apresentadosanteriormente.Consideremosentão a seguinte aplicação:(Ç,1]) E fib t-i (x,y) E fiedefinida porx(Ç,1]) 3LcPa(ç,1])X a l15(50)
3L Pa([:.,11)Y:a lY([:.,11) (51)As funções (50) e (51) são biunívocM entre o triângulo ne e o triânguloretângulo nb.Definindo([:.1,111) (0,0)([:.2,112) (1,0)3{Xa.e0.sesea a be{ Pa([:.b,11b)Ya Ya'e0,sesea ba beX([:.b,11b) Pa([:.b,11b)Xa Y([:.b.11b) L3([:.3,113) (0,1)(52)ba lMostramos a transformação geometricamentena Figura 6, onde y y( ) (x.Y) e (y) ([:.,11).(0,1) .(X3.y3)y.o.b-1(XI.Y (0,0)(1,0)(x y;)Figure 6: Thansformação IsoparamétricaPara verificar a existência da funçãoy-1 : ne H- Obusaremos o Teorema da F\1nção Inversa. A função Paé diferenciá.vel e portantoy( ) é diferenciável. Assim podemos calcular o Jacobiano da transformaçãoisoparamétrica, dado por 8zJ det 8zt 08F. 81/Como y : nb H- ne é bijetora e de clMse C1 e o Jacobiano é positivo, então, oTeorema da F\1nção Inversa garante a existência da função y-1 : ne H- Ob,16
com y-1 de classe CI. Assim temos um mapeamento entre os elementos finitosnb e ne.A função y( ) pode ser representada da seguinte forma:x( , '7) (x -x )Ç (x -x )'7y(Ç, '7) (y -yi)Ç (y -yi)'7 xi yi(53)Cálculo do JacobianoDe y( ) definido em (53), podemos calcular o Jacobiano da transformaçãoentre nb e ne, ou seja 8zJ JA função det 8 8'1(X -xi)(Y de interpolação[ (X2e8z 8(" det(e-XIe -eY2-yi)-(x YI))-x )(Y2( ex3 -XI( e -eY3{ Pa( ) PI( ) 1 -ç P2( ) çYI))]-y ) Pa( ) pode ser representadaformae(54)em nb na seguinte-1](55) P3( ) 1]Daí, temos que a derivada da função Pa( ) é dada por -18Ç 8 , -181] 18Ç, o81} o8Ç, 18'7(56)Assim, o gradiente da função Pa( ), é dado por{V' pa( ) V' PI( ) (-1,-1)V' P2(Ç) (1, O)V' P3( ) (0, 1)(57)Temos que :ne(x,y)-tH-nb (x,y) (Ç,l])A função (Y) pode ser dada comoÇ(x, y) J [ (Y -Yi)x (x -x3)Y (x3Yi -x Y ) ](58)1J(x, Y) J [ (Yi -Y )x (x -xi)y (xiY2 -x Yi)(59)17]
Sendo (y) uma função linear que satisfazf.(Xl,Yl) O17(Xl,Yl) Of.(X2,Y2) 117(X2,Y2) Of.(X3,Y3) O1}(X3,Y3) 1Com as matrizes globais M, C, K, e o vetor força global G montados apartir das matrizes locais e do vetor força local é possível, agora, escrever (36)numa forma matricial dada por:Md(t) cd(t) Kd(t) G(60)onde a matriz de massa global M é simétrica e a matriz de amortecimento globalC é anti-simétrica. Sendo d [dl , d2, ., dm Jt O vetor incógnita.8MétododeN ewmarkMatematicamente, a equação (60) representa um sistema de equações diferenciais ordinárias de segunda ordem e, em particular , a solução das equações nãosão simples já que as matrizes dependem de y em cada instante de tempo.E os procedimentos usualmente propostos para a solução de sistemas gerais deequações diferenciais podem ser muito caros se a ordem das matrizes for grande.O Método de Newmark foi escolhido para resolver o sistema de equaçõesdiferenciais de segunda ordem representado pela equação (60 ) .Este método éum Método de Diferenças Finitas que foi apresentado por Newmark em 1959, ea partir de então vem sendo largamente utilizado para a resolução de problemas.dinâmicos.O método de Newmark nos dá a solução no final de um passo do tempoexpressada por uma série de Taylor aproximado por:d(t 6.t) d(t) 6.td(t) [ (1-2f.)d(t) 2f.d(t 6.t) ]d(t 6.t) d(t) 6.t [ (1 -15)d(t) 15d(t 6.t) ](61)(62)onde f. e 15são parâmetros livres que determinam a estabilidade e precisão características do algoritmo considerado. A especificação dessas constantes nosfornece uma variedade de métodos conhecidos como Família de Métodos deNewmark.O problema de valor inicial para (23) consiste em encontrar um deslocamentod(t) tal que,Md(t) cd(t) Kd(t) GIV d(O) do(63)d(o) do{18
onde d(O) e d(o) são os deslocamentos e as velocidades iniciais dos nós da cordaconhecidos. As equações (63) formam o que chamaremos de Problema (IV).Do Problema (III), (23), temos ainda as condições de fronteira dadas pelasegunda equação de (23) que devem ser respeitadas, ou seja,d(t) Opara y O e y 1(64)Como d(t) corresponde a velocidade nodal e d(t) a aceleração nodal, faremosa seguinte associação d(t) v(t) e d(t) a(t). Assim de (60), (61) e (62) temos,para cada instante de tempo tn nÂt, o seguinte sistema iterativo:M an l a'IJn l K dn l GÂt2dn l dn Âtvn 2[(1-21:.)an 2Çan JVn l Vn Ât[(1 -d)an(65) dan l]onde an, Vn e dn representam d(tn), d(tn) e d(tn) respectivamente.De (65) temos três equações para determinarmos três incógnitas an l, 'IJn ledn l.Neste trabalho, escolhemos trabalhar com I:. 0 e d o que correspondea um Método de Diferença Central. Assim (65) pode ser reescrita como:Man l aVn l Kdn l G(66)Ât2dn l dn Âtvn 2an(67)Ât'IJn l Vn (an(68) an l)ImplementaçãoExistem várias possibilidades para a implementaçãoescolhemos a seguinte forma.deste método, assim-Dados IniciaisPara tn 0, as matrizes de massa global M, de amortecimento global a erigidez global K, e o vetor global G são montados como foi apresentado anteriormente utilizando o Método dos Elementos Finitos.De (63) temos que o deslocamento inicial e a velocidade inicial Vo sãodados. A partir daí é possível calcular a aceleração inicial ao diretamente de(66), ou seja:Mao G- avo -Kdo(69)resolvendo o sistema linear (69) nós obtemos a aceleração inicial ao.-Demais instantes de tempoPara cada tn, as matrizes globais M, a e K, e o vetor global G são montadoscomo foi apresentado utilizando o Método dos Elementos Finitos.19
SubstituindoMan l O(67) e (68) em (66) temos(Vn 2(anÂt an l)) K (dn Âtvn TanÂt2) G(70)Separando em (70) no lado esquerdo da equação o que for relativo ao tempotn l e do lado direito o que for relativo ao tempo tn, obtemosMan l(M o2an lÂt 02ô.t) an l G-O G - (Vn 2anÂt) -K (dn tvn)(71)(Vn 2an ) -K (dn Ô.tvn 2an )(72)ô.t 2anÂt2ô.t2Definimos os preditores como sendo:-ô.t2dn l dn Âtvn 2an(73)iin l Vn an(74)Logo as equações (67) e (68) podem ser reescritas como:dn l dn l(75)Vn l iin l an l(76)Substituindo (73) e (74) em (72) temos que a relação recursiva para determinar an l é dada por(M 2 ) an lÂt Gn l -oiin l-Kdn l-(77)As equações (75) e (76) serão utilizadas para calcularmos dn l e Vn l, respectivamente.Logo, para cada instante de tempo tn nÂt com 0 . : tn T obtem-se dn,com n E 1N.Da equação (27) temos quemvh(y, t) 2::: di(t) pi(Y),i l PiE VmCom os deslocamentos nodais dn l calculados encontramos Vh(y,t) dadopela relação (27). Assim temos a solução do Problema (111).Mas, vimos anteriormente que resolver o Problema (III) implica em encontraruma solução aproximada para o Problema (II). E como a transformação A, dadaem (14), é um isomorfismo podemos afirmar que resolver o Problema (II) éequivalente a resolver o Problema (1).20
9SoluçãoGlobalPara a obteção da solução global v(y, t), para o problema cilíndrico definidoem (13) Medeiros e FerreI introduziram uma viscosidade.5v'(y, t), onde .5é umaconstante positiva fixada, ou seja, que a força deveser proporcional a velocidadedo deslocamento. Assim, o problema penalizado é dado pelo seguinte problemacilíndrico:{V Lv(y, t) .5v'(y, t) g(y,t),(y, t) E Qv(y, t) O,y EE8vv(y,O) vo(Y),at(y,O) Vl(Y),(78)yEnonde Vo e Vl são conhecidose Lv(y, t) é o operador definido em (16)Lv(y,t)82V -82v82v8Yi 8 Yj Ci(Y,t)- 8t8 Yi8 t 2 -ã(t)6.v bij(Y,t)8v ei(Y,t)-8 Yiondeã(t) b ( K'(t) ) YiYjCi(Y, t) -2ei(Y,t) rl2bij(Y,t)( a(t) b(t) L I\7VI2dy)( K'(t) ) Yi( K'(t) ) 2 K"(t) 1-K(t)JYiA introdução da viscosidade .5v'(y, t) se faz necessária para o estudo dasolução global. Em [4] Medeiros e FerreI fizeram toda a sua análise do modeloutilizando essaviscosidade.Desde que u(x,t) v(y,t), comxy K(t)então: - 8t(-K(t) .x1) 8Xi 8Xi 8tpor conseqiiênciaLA (U x,t) ar[( K(t) x, ) .21 8t] i( x, t )(79)
Assim, o seguinte problema não-cilíndricoé obtido, equivalente ao problema(78):{ Lu(x,Au(x,t)u(x,o)t) 8[() âu(x,t)KI(t)K(t)Xi âU(X,t)ât O] f(x, t)A(x, t) E Q(x,t)âu uo(x), "8"i(x, 0) Ul(X)E Ê(80)x E noPortanto, encontrando a solução global de (78), obtém-se uma solução globalpara o problema misto (80) associado M vibrações de uma membrana elásticacom extremidades móveis, devido a equivalência dos problemas mostrado naseção 2.10SimulaçõesN uméricasExemplos numéricos serão mostrados nesta seção para ilustrar algumas características do modelo para pequenas vibrações transversais de membranas elásticascom fronteira livre.Vimos que resolver o Problema III, (23), ou seja, encontrar vh(y, t) implicaem encontrar uma solução aproximada do Problema (II), (21), e assim resolver oProblema (I), (13), ou seja, encontrar u(y, t) já que o Problema (I) é equivalenteao Problema (II).Para os exemplos numéricos consideraremos uma membrana elástica subdividida em 73 nós e 128 elementos finitos e o passo do tempo utilizado serádt 0,0025.10.1Exemplo1Consideraremos a posição inicial da membrana e a velocidade inicial dadas por:vh(y,o) 0(81) (l-x2-y2).!(82)ât7rO crescimento do raio da membrana, apresentado no gráfico 7, é dado por:K(t)t 35 1(83)Assim como para as cordas elásticas a força externa aplicada na membranaelástica é normalmente considerada nula, mas primeiramente nós construiremos um exemplo numérico onde a força externa não é nula, com o intuito deconfirmarmos que a solução está sendo obtida corretamente.22
tadoProblemaIIIcomosendo:1vh(y,t) (1-X2-y2)2sen(11"t)(84)11"ObservemosParaquea �sticasfísicasTO- 5eminiciaissubstituirmos(84)daem(81)e (82).(23).membrana:k- ovimentoeaproximadaIIIestádamem-damem-nula.oé emoscorretamente,aO, �central
O o0040020-0.02-0.04-O.-O. 0FigureApesaratravésesperadoEmIado modelodo gráficoestar sendo resolv do corretame8, quandopara a membrana[4], Medeirostomamosem vibre FerreIa força-nãofizeramexterte, o q e foi comprovadoa nuloacontece comoda a análiseo modvibrações transversaisde membranasel ticas com fr nteiraa penalizaçãoapresentadana seção 9. Assim, enco tramoexemplo 1 considerandoa força externnula quando tomestá. apresentadono gráficoamortecimentoostraO gráfico9.10 para pequenasivre considerandoa solução para oos ó O, 5, o que10.o. o.m0.040.020-0.02-0.04-O. 0Figure10.2ExemploConsideraremos810: v h (0,0, t) qu10do a força ex erna éuIa2a posiçãoinicialda mebranae a vel cidadinicialdadas por:Ivh(t/,O) (1X2-y2)2(85)11"
oat( 86 )O crescimento do raio da membrana, apresentados no gráfico 11, é dado por:K(t)t.7tAt.stAt.3t t.t /t ,o 2- e-o,OSt(87) .///K(I)/101520Figure 11: Crescimento do raio da membranaO gráfico 12 apresenta o movimento do nó central da membrana quando aforça externa é nula, !5 1 e as característicasfísicasda membrana consideradascomo:TOk- 10 e - 1mmOt2 I0.O 0040020.002.004.o.o.Q10tOFigure 12: u(o, 0, t) quanto tomamos !5 1Apresentamos na Figura 13 as posiçõesda membrana para os tempost o;0,25; 0,5; 0,75; 1; 1,2525
l 001o":Oo1 Cl":Lo,.--0:;;:-r s--: ?',-Ol 001o. O o.o.,.:d;s---r s--:; 1-Oo.01.1."-45 'l1-0,25 1.-ç 'n- :("n- :("'-G,5Cl":001l1-.7' 00": Oo1L "' .,1 ?'.:d;s---r s-- -1 :SCJ"'"'AI-I--1. --r s-- -1-.:SCJ"'"'A1I -1,25Figure 13: u(y, t) quanto tomamos 6 110.3Exemplo3Consideraremosa mesma posição inicial e a masma velocidade inicial da membrana consideradasno exemplo 2 dadas em (85) e (86):vh(y,O) (1- X2 -y2)2111" OâtO crescimentodo raio da membrana, apresentadosno gráfico 14, é dado por:K(t) 1,1 -0,1cos (Kt)(88)Nos exemplos anteriores tomamos K(t) crescente. Verificaremos através docorrente exemplo o que ocorre quando tomamos uma função K(t) que não ésempre crescente.Observamos, através do gráfico 14, quando t varia de O a 5 instantes detempo a membrana cresceaté atigir o seu tamanho máximo em t 5 quando26
1.2/ cresceéemnula,6 tingir150,dootseueasvariadetamanhoo ranadatempoaquandomembranaa1.a forçaconsideradascomo:.:-1!. 8me! vido,deovibraçãoquedanãoclaramentejáem16a blemaI.27K(t)quenãoésemprexy
020,1- '008008004-,0.5/ ,-'--- :;---G:S u --r--:1--- .o:s5., ;0.5 o--0----;0 U---"", .,'-0.012 --- ,008O 0040020'-- , .-"-"1., / .r.C.:;--- --- ( 5-0 -T-Õ:s---r- 0--/;0.5 J',-,0 1'2008O l0/- /á -:,/'--"- ;::----,-,.(. --- .r-.a:s--- ./;/-,/.\/.5/,;-:o- u --r--/.1, -10 .0.'"67,/'r 5-0---. ""- 1./ .'- Figure 16: u(Y. t) e sua projeção no plano xyReferences[1] CARR1ER. C. E.j On the vibration problem of elastic string .Q. J. Appl.Math. 151-165 (1953)[2] HUGLES. T. J. R.j The Finite Element Method Linear Static andDynamicFinite Element Analynis. Prentice Hall (1987)[3] MEDEIROS. L. A.j FERREL. J. L.j Kirchhoff-Carrierelastic string in noncilindrical domain
extensão do modelo de Kirchhoff-Carrier, uma vez que leva em conta a troca de tamanho da membrana durante a vibração. O modelo de Medeiros e Ferre] não possui solução analítica, então usamos o Método dos Elementos Finitos para a resolução do modelo no espaço e um Método de Diferenças Central, o Método de Newmark, para a .
