Practica 2 Oscilaciones con 1 grado de libertad Catedra de Mecanica. 16 de Noviembre de 2004. Practicas de Mecanica Computacional 1. Oscilaciones forzadas con amortiguamiento > restart; > with(plots): Warning, the name changecoords has been redefined

> with(plottools): Warning, the name arrow has been redefined

Se define la ecuacion diferencial y se asigna a oscilforc: > oscilforc:=m*diff(x(t),t$2)+c*diff(x(t),t)+k*x(t)=q*sin(Omega *t);  d2  d  oscilforc := m  2 x(t ) + c  x(t ) + k x(t ) = q sin(Ω t )  dt   dt  Con la instruccion dsolve se obtiene la solucion general de la ecuacion diferencial: > dsolve(oscilforc); x(t ) = e

  (c −  − 

 c2 − 4 k m ) t   2m 

_C2 + e

  (c +  − 

 c2 − 4 k m ) t   2m 

_C1 +

q ((k − Ω 2 m) sin(Ω t ) − Ω cos(Ω t ) c ) Ω 4 m 2 + (c 2 − 2 k m ) Ω 2 + k 2

Con la instruccion assume se impone la condicion de que el amortiguamiento sea subcritico: > assume(c**2 solgen:=dsolve(oscilforc); solgen := x(t ) = e +e

 c~ t   −   2 m~ 

 c~ t   −   2 m~ 

 −c~2 + 4 k~ m~  sin  2 m~

 −c~2 + 4 k~ m~  cos  2 m~

t   _C2 

t  ((−k~ + Ω 2 m~ ) sin(Ω t ) + c~ Ω cos(Ω t )) q  _C1 −  Ω 4 m~ 2 + (−2 k~ m~ + c~2 ) Ω 2 + k~2

Observese como las variables m, c y k, a las que se les ha impuesto una restriccion mediante assume en lo sucesivo se denotan seguidas de una tilde: m~, c~ y k~

En la siguiente instruccion se define la funciónxtran aplicando la instruccion unapply, con el argumento t para indicar que es funcion del tiempo, al lado derecho de la igualdad (operador rhs) que se ha asignado anteriormente a solgen:

> xtran:= unapply(rhs(solgen),t): A continuacion se define la funcion correspondiente a la solucion de regimen permanente. Para ello, con op(3,xtran(t)), se selecciona el tercer sumando de xtran(t), y con unapply(*,t) se establece que la funcion depende del tiempo: > xperm:=unapply(op(3,xtran(t)),t): --------------------------------------------------------Aplicacion al ejercicio de practicas 15: Los datos del oscilador son: m=400, c=368.4136, k=100000 Los datos de la fuerza externa: q=3158.2734, Omega=4*Pi y los dos conjuntos de condiciones iniciales: a) x0=0.08, v0=0 b) x0=0, v0=1 --------------------------------------------------------Se vuelve a resolver la ecuacion diferencial pero considerando ahora en dsolve las correspondientes condiciones iniciales. En solejA se asigna la solucion que se obtiene considerando el conjunto (a) de condiciones iniciales: > solejA:=dsolve({oscilforc,x(0)=0.08,D(x)(0)=0},x(t)): y se define xA como la funcion del tiempo correspondiente a la solucion solejA empleando otra vez la instrucción unapply: > xA:= unapply(rhs(solejA),t): Analogamente, en solejB se asigna la solucion que resulta de considerar el conjunto (b) de condiciones iniciales: > solejB:=dsolve({oscilforc,x(0)=0.0,D(x)(0)=1.0},x(t)): y en xB se define la funcion del tiempo correspondiente a esta solucion: > xB:= unapply(rhs(solejB),t): A continuacion se particularizan las funciones xA y xB para los valores numericos del enunciado. Para ello se sustituyen, mediante el operador subs, los valores de m, c, k, q y Omega en xA(t) y con el resultado obtenido se define una nueva funcion del tiempo que llamamos xanum y xbnum, respectivamente: > xanum:=unapply(subs(m=400,c=368.4136,k=100000,q=3158.2734,Ome ga=4*Pi,xA(t)),t): > xbnum:=unapply(subs(m=400,c=368.4136,k=100000,q=3158.2734,Ome ga=4*Pi,xB(t)),t): Finalmente, con la instruccion plot se dibujan las ecuaciones horarias xanum(t) y xbnum(t), pudiendo comprobar que si bien los regimenes transitorios correspondientes a ambas soluciones son diferentes, el regimen permanente no depende de las condiciones iniciales. > plot({xanum(t),xbnum(t)},t=0..8);

0.15

0.1

0.05

0

1

2

3

4

5

6

7

8

t –0.05

–0.1

> a:=animatecurve(xanum(t),t=0..8,numpoints=200,frames=50,color =blue,thickness=3): > b:=animatecurve(xbnum(t),t=0..8,numpoints=200,frames=50,color =red,thickness=3): > display(a,b): > dibu:=proc(s) > r1:=rectangle([-1,-10],[0,10],color=black): > r2:=rectangle([0,-1],[s,1],color=blue): > d1:=disk([s,0],3,color=red): > display([r1,r2,d1]): > end: Warning, ‘r1‘ is implicitly declared local to procedure ‘dibu‘ Warning, ‘r2‘ is implicitly declared local to procedure ‘dibu‘ Warning, ‘d1‘ is implicitly declared local to procedure ‘dibu‘

> vals:=seq(xanum(i*2/50),i=0..50): > fotos:=map(u->dibu(u*500),[vals]): > display(fotos,insequence=true,scaling=constrained);

10 –60

–40

–20

0

20

40

60

–10

Para dibujar las ecuaciones horarias no es necesario definir una funcion mediante unapply. Lo que si hay que hacer es definir a priori los valores numericos: > mp:=400:cp:=368.4136:kp:=100000:fp:=3158.2734:omp:=4*Pi:x0p:= 0.08:v0p:=0: > solpru:=dsolve({mp*diff(x(t),t$2)+cp*diff(x(t),t)+kp*x(t)=fp* sin(omp*t),x(0)=x0p,D(x)(0)=v0p},x(t)): > plot(rhs(solpru),t=0..8);

0.15

0.1

0.05

0

1

2

3

4

5

6

7

8

t –0.05

–0.1

2. Frecuencia y Amplitud de Resonancia Para hacer el calculo de la frecuencia y de la amplitud de resonancia expresamos la solución de regimen permanente como el producto de la amplitud A por una funcion coseno con un cierto desfase a. Dicha solucion la asignamos a la funcion del tiempo xperm2: > xperm2:= t -> A*cos(Omega*t+alpha); xperm2 := t → A cos(Ω t + α ) Para obtener la expresion de A y a identificamos los valores de las funciones xperm (definida mas arriba) y xperm2 en los instantes de tiempo correspondientes a t=0 yt=p/2W: > ec1:=xperm(0)=xperm2(0); ec1 := −

c~ Ω q Ω 4 m~ 2 + (−2 k~ m~ + c~2 ) Ω 2 + k~2

= A cos(α )

> ec2:=xperm(Pi/(2*Omega))=xperm2(Pi/(2*Omega)); ec2 := −

(−k~ + Ω 2 m~ ) q Ω 4 m~ 2 + (−2 k~ m~ + c~2 ) Ω 2 + k~2

= −A sin(α )

Resultando (a la amplitud de regimen permanente la denominaremos Aux1 para manipularla posteriormente): > alpha:=arctan(lhs(ec2)/lhs(ec1));  −k~ + Ω 2 m~   α := arctan  c~ Ω  > Aux1:=simplify(solve(ec1,A)); k~2 − 2 k~ Ω 2 m~ + Ω 2 c~2 + Ω 4 m~ 2

q c~ Ω Aux1 := −

Ω 2 c~2

k~2 − 2 k~ Ω 2 m~ + Ω 2 c~2 + Ω 4 m~ 2

La solucion obtenida de la amplitud Aux1 la simplificamos usando la instruccion simplify con la opcion symbolic, para obtener una expresion similar a la de la ecuacion (3.22) del libro de apuntes: > Aux2:=simplify(Aux1,symbolic); q

Aux2 := −

k~2 − 2 k~ Ω 2 m~ + Ω 2 c~2 + Ω 4 m~ 2

La expresion de la amplitud de resonancia la obtendremos expresando la amplitud de regimen permanente Aux2 en funcion de la frecuencia W de la fuerza externa: > Amp:=unapply(Aux2,Omega); q Amp := Ω → − k~2 − 2 k~ Ω 2 m~ + Ω 2 c~2 + Ω 4 m~ 2 y calculando su maximo. Para ello, derivamos respecto de W empleando la instruccion diff(Amp(Omega),Omega), igualando a cero, > wres_eq:=diff(Amp(Omega),Omega)=0: La ecuacion anterior se resuelve con la instruccion solve: > wres_sol:=solve(wres_eq,Omega); wres_sol := 0,

4 k~ m~ − 2 c~2 ,− 2 m~

4 k~ m~ − 2 c~2 2 m~

De las tres soluciones obtenidas, la que tiene sentido fisico es la segunda: > wres:=wres_sol[2]; wres :=

4 k~ m~ − 2 c~2 2 m~

A continuacion podemos hacer una manipulacion puramente "cosmetica" para conseguir la expresion (3.25) del libro de apuentes: > expand(wres**2)**(1/2); k~ c~2 − m~ 2 m~ 2 La amplitud de resonancia se obtiene sustituyendo la frecuencia de resonancia wres en la

La amplitud de resonancia se obtiene sustituyendo la frecuencia de resonancia wres en la expresion de la funcion Amp: > A_reson:=Amp(wres): Esta expresion la simplificamos: > A_reson:=simplify(sqrt(factor(Amp(wres)**2)),symbolic); 2 q m~ A_reson := c~ −c~2 + 4 k~ m~ Con objeto de representar la amplitud en funcion de la frecuencia y del amortiguamiento se adimensionaliza la funcion definida anteriormente Amp(Omega) en terminos de la fraccion de amortiguamiento critico x y de la frecuencia natural w. Para ello se sustituye con la instruccion subs, el parametro c por sqrt(4*k*m*xi), > A_aux1:=subs(c=sqrt(4*k*m*xi),Amp(Omega)); q A_aux1 := − k~2 − 2 k~ Ω 2 m~ + 4 Ω 2 k~ m~ ξ + Ω 4 m~ 2 > A_aux2:=sqrt(A_aux1**2*(k/q)**2): se sustituye el parametro k por omega**2*m, > A_aux3:=simplify(subs(k=omega**2*m,A_aux2)); A_aux3 :=

ω4 ω4 − 2 ω2 Ω2 + 4 Ω2 ω2 ξ + Ω4

y finalmente se sustituye b=W/w mediante la instruccion subs(Omega=beta*omega,A_aux3) > A_aux4:=simplify(subs(Omega=beta*omega,A_aux3)); A_aux4 :=

1 1 − 2 β2 + 4 β2 ξ + β4

Para dibujar las curvas se define la funcion f(x,b) a partir de la expresion asignada a A_aux4: > f:=unapply(A_aux4,beta); f := β →

1 1 − 2 β2 + 4 β2 ξ + β4

y se consideran los casos correspondientes a las fracciones de amortiguamiento critico x=0.0, x=0.05, x=0.20, x=0.50, x=0.70 y x=1.0 > f_00:=subs(xi=0.0,f(beta)): > f_05:=subs(xi=0.05,f(beta)): > f_20:=subs(xi=0.20,f(beta)): > f_50:=subs(xi=0.50,f(beta)): > f_70:=subs(xi=0.70,f(beta)): > f_100:=subs(xi=1.0,f(beta)): Las seis funciones se dibujan en un mismo grafico, agrupandolas entre llaves, con la instruccion plot > plot({f_00(beta),f_05(beta),f_20(beta),f_50(beta),f_70(beta),

> plot({f_00(beta),f_05(beta),f_20(beta),f_50(beta),f_70(beta), f_100(beta)},beta=0..3,y=0..10);

10

8

6 y 4

2

0

0.5

1

1.5

2

2.5

3

beta

3. Respuesta a cargas arbitrarias definidas "a trozos" Tambien es muy sencillo imponer cargas arbitrarias definidas "a trozos". Consideremos la accion de tres cargas triangulares: trian1: actua durante un periodo de tiempo que es aproximadamente la mitad del periodo propio del oscilador trian2: actua durante un periodo de tiempo aproximadamente cinco veces el periodo propio trian3: actua durante un periodo de tiempo veinte veces superior al del oscilador > Tprop:=evalf(2*Pi/sqrt(kp/mp)); Tprop := 0.3973835307 > triang1:=PIECEWISE([1000*t, t triang2:=PIECEWISE([100*t, t soltrian1:=dsolve({mp*diff(x(t),t$2)+cp*diff(x(t),t)+kp*x(t)= triang1,x(0)=0.0,D(x)(0)=0},x(t)): > soltrian2:=dsolve({mp*diff(x(t),t$2)+cp*diff(x(t),t)+kp*x(t)= triang2,x(0)=0.0,D(x)(0)=0},x(t)): > soltrian3:=dsolve({mp*diff(x(t),t$2)+cp*diff(x(t),t)+kp*x(t)= triang3,x(0)=0.0,D(x)(0)=0},x(t)): > plot([rhs(soltrian1),rhs(soltrian2),rhs(soltrian3)],t=0..10);

0.001

0.0005

0

2

4

6

8

10

t –0.0005

–0.001

> > > > >

4. Ejercicio propuesto Para el oscilador cuyas constantes se han empleado a lo largo de la practica (m=400 kg, c=368.4136, k=100000) obtener y representar la respuesta para tres cargas escalon de valor F=1, que comienzan a actuar en t=0, y que se aplican durante 0.2, 10 y 15 segundos respectivamente.