45

• Capítulo 3- Sistemas Vibratórios Excitados por Fonte de Potência Limitada

Pesquisando-se a vasta bibliografia disponível, consta-se, contudo que é usual

desconsiderar-se, no modelo matemático, a influência do movimento do próprio sistema em sua excitação. Todavia, em muitos casos, esta simplificação, no modelo matemático, não é razoável, devendo-se levar em conta que a excitação ou sua fonte, é influenciada pela própria resposta do sistema.

Este fato prejudica a formulação dos modelos matemáticos da Teoria de

Vibrações tradicional, necessitando-se de estabelecer uma formulação mais realística que leve em conta a interação entre as variáveis de controle, e as de estado, da excitação, com as de estado do sistema físico estrutural.

Tem-se assim, um sistema vibratório não-ideal, ou um sistema com fonte de

excitação, não-ideal. O sistema ideal é o tradicional onde não existe este fenômeno. Note que o sistema não-ideal, tem sempre um grau de liberdade a mais do que seu ideal, correspondente.

Entenda-se, que a máquina não-ideal é uma conceituação que depende

fundamentalmente da estrutura que a suporta. Um exemplo clássico de um sistema não-ideal é o de uma estrutura flexível (por

exemplo uma viga em balanço)sobre o qual um motor elétrico de corrente contínua desbalanceado, com fonte de potência limitado, está montado.

Nestes casos, mencionados, o movimento do sistema devido a sua própria flexibilidade, afeta o desempenho da máquina. Notam-se, características ,importantes, deste tipo de modelo matemático de sistema dinâmico, provenientes da interação da estrutura com a fonte de excitação, tais como:

• Variações bruscas ( saltos ou “jumps” ) da amplitude de deslocamento da estrutura e da freqüência de excitação, no caso particular em que considera-se as condições estacionárias do movimento;

• Descontinuidade da curva “amplitude versus freqüência “;

• Influência do perfil da curva “amplitude versus freqüência “ quando o operador

altera a velocidade do motor elétrico, em acréscimos e/ou decréscimos;

• Dependência destes efeitos com as características eletromecânicas do motor.

Portanto, percebe-se que os sistemas dinâmicos, modelados matematicamente, como sendo não-ideais, possui um grau de liberdade(ou mais dependendo do número de motores presentes ,no sistema) ,superior ao sistema ideal correspondente.

46

O fenômeno, em tela, é conhecido, na literatura como efeito Sommerfeld, em homenagem ao primeiro pesquisador, a observar estes fatos, experimentalmente(Sommerfeld,1903). Recentemente este problema foi analisado por nós, constatando que as não-linearidades quadráticas e cubicas, deste problema, tinham a mesma ordem de grandeza( De Mattos et al., 1997). Tal fato foi importante na modelagem segura deste problema de vibrações não-lineares. Recentemente, demonstrou-se haver boa concordância entre os resultados numéricos e experimentais(De Mattos et al., 1999). Estes 2 trabalhos, deram a este problema, avanço significativo, na literatura. Nota-se que as equações governantes do movimento, representativa do sistema dinâmico não-ideal de vibrações, podem se esquematizadas, modelados matematicamente, da seguinte forma:

1-“ primeiro conjunto de equações ”:

equações dinâmicas da “estrutura vibratória ”em x + termos interação ( termos em x e

em ϕ ) = equações da excitação externa.

2-“ segundo conjunto de equações”:

equações de rotação do eixo do motor ( termos em 2

2

d

td

ϕ ) + equações características do

motor ( termos em L ( d

d t

ϕ ) )+ outros termos de interação( em x e em ϕ ) = 0.

onde

• d

d t

ϕ é a velocidade angular do eixo do rotor e,

• 2

2

d

td

ϕ a aceleração correspondente ,onde supôs-se a existência de um único

motor. Note o fato importante de que as linearidades ou não-linearidades vem da

modelagem matemática da estrutura, no caso do problema ideal; a presença de por exemplo, uma única fonte de energia, em geral contribui com a ação de não-linearidades do tipo quadráticas, presente no termo de energia cinética, devido o desbalanceamento do motor de corrente contínua (Nayfeh e Mook, 1979).

Ressaltam-se que as vibrações oriundas dos movimentos de aceleração (partida) e

desaceleração (desligamento) de motores elétricos são representados pela função L, ou seja ela representa o torque desenvolvido pelo motor elétrico e a curva característica do motor

47

relaciona o torque L com a velocidade de rotação angular. Cada ponto ,desta curva,

considerada, pode ser obtida da seguinte maneira: Mantém-se a velocidade angular d

d t

ϕ

constante, enquanto se determina o torque no eixo do motor. Como cada ponto gira a uma velocidade angular constante, essa curva é chamada de característica estática, determinando seu funcionamento em estado estacionário. Na realidade, todo motor possui uma família destas curvas, em discussão, e todas elas representam a mesma relação de grandezas, mas a cada uma corresponde um certo valor de Regulagem.

Imagine-se, por exemplo, que um motor possua uma família de curvas, as quais são diferenciadas variando-se a tensão elétrica U que é aplicada aos seus terminais. Cada posição de controle de regulagem com o qual se aumenta ou se diminui a tensão elétrica, produz uma correspondente curva característica, relacionadas, cada uma delas, a um valor de potência constante, mas diferentes entre si.

As Figura 1 e 2, ilustram duas possibilidades destes gráficos, no plano torque

versus freqüência de rotação.

• A primeira delas sendo do tipo exponencial(mais realística dt

dE

eEL

ϕ2

1

−

= ) e;

• A segunda sendo do tipo linear(dt

dbaL

ϕ+= ) para cada valor da tensão U.

U será o parâmetro de controle do problema. A seguir, mostra-se a dedução das

Retas características; as curvas características, tem determinação experimental.

Figura 1-Curvas características Figura 2-Retas características A seguir, explicita-se a obtenção das retas características, desde que as curvas

características tem determinação exclusivamente, experimental. Admite-se motores de corrente contínua esquematizado, através da Figura 3, cujas

equações de controle representativas de seus circuitos elétricos, são dadas pela equação, clássica:

U-ea=Raia+Laia

Onde:

48

• U é tensão elétrica aplicada aos terminais do motor;

• ea, a fôrça contra-eletromotriz;

• Ra a resistência elétrica do motor;

• ia sua corrente elétrica e;

• La sua indutância elétrica.

Figura 3: Modelo de Armadura do Motor

Tem-se:

ea = ke &ψ e TE = kmia

com (ke, km, TE, &ψ ): constante de tensão do motor, constante do torque, torque eletrotromecânico e velocidade angular.

A lei que vincula variação do torque, torque aplicado, tensão elétrica no motor e a

velocidade angular é da forma:

Uke &ψ =(Ra/km)TE+ (La / km) TE

Em Regime Estacionário, a variação temporal do torque TE é zero(ouLa é

desprezível), logo tem-se:

TE=U(km/Ra)-ke(km/Ra) &ψ

Nota-se, então se a tensão U for mantida constante, a relação existente entre o

torque e a velocidade angular é do tipo linear, portanto as curvas características deste motor, são retas interceptando os dois eixos, como mostrado na Figura 1.

Pode-se modelar matematicamente os torques, em questão, por:

TE=Tm+I &ψ

com

49

• TE sendo o torque eletro-mecânico e;

• Tm o torque mecânico e I o momento de inércia do eixo do motor.

O torque mecânico Tm pode ser modelado como sendo da forma (Hindmarsh, 1981):

Tm=k1+k2 &ψ +k3( &ψ )2+T(t)

Onde:

• O primeiro termo está ligado ao atrito constante do tipo Coulomb devido ao atrito dos mancais;

• O segundo termo é o atrito viscoso associado a não turbulências e;

• O terceiro termo é devido a turbulência, i.e., a ação da propulsão ou dos

ventiladores refrescantes e o último é o carregamento da estrutura que suporta o motor.

Por conveniência e por facilidade de comparações com a literatura, adota-se daqui para

frente, a notação:

{ TE=L, &ψ =x4, M(x4)=L(x4)-H(x4),

onde

H=g1 x4 e ψϕ = }

Nos próximos capítulos, definem-se os modelos matemáticos, que serão objetos de

estudo deste trabalho.

Finalmente, mencionam-se ,algum conceitos ,que valem tanto para sistemas dinâmicos ideais e não-ideais, na região de ressonância, isto é, definida por:

d

d t

ϕ - ω = Ο ( ε1 )

onde: ε1 é um pequeno parâmetro e, ω é a freqüência natural do sistema dinâmico vibratório. Supondo-se que o sistema dinâmico, em questão, partisse do repouso, a velocidade

angular d

d t

ϕ do eixo do rotor cresceria até atingir a região de ressonância.

Assim, dependendo das condições iniciais e dos parâmetros físicos, tais como;

50

• Massa desbalanceada;

• Massa do motor;

• Momento de inércia do rotor;

• Excentricidade da massa desbalanceada;

• Constantes de torque desenvolvido pelo motor:

A tensão aplicada nos terminais do motor, d

d t

ϕ continuaria crescendo além da

região de ressonância(Fenômeno da Passagem pela Ressonância)ou, Permaneceria próximo a freqüência natural ω do sistema(Fenômeno da Captura pela Ressonância). O tempo de passagem pela ressonância, dependerá, também das condições iniciais do sistema.

O primeiro autor a sistematizar estas propriedades dos modelos matemáticos , não-ideais, foi( Kononenko, 1969). Atualmente, um Completo e compreensível revisão deste assunto, considerando o período de 1904-1999, pode ser encontrado em (Balthazar et al., 1999). O fenômeno de Salto e o aumento de potência exigido pela fonte de energia, operando na região da ressonância são manifestações do fenômeno, conhecido como SOMMERFELD.

Nota-se também uma característica importante existente em sistemas não-ideais: a maioria das máquinas rotativas mostram características de decréscimo de torque quando é operada com acréscimos de velocidades de rotação.

Finalmente, com o objetivo, de esclarecimento final, exibe-se, na Figura 3 , um

exemplo típico da relação existente entre a velocidade de rotação o torque atuante(Yamakaswa e Murami, 1989).

Na Figura 4, notam-se as características do torque atribuído à máquina, mostrado

por TS e a relação existente entre a velocidade de rotação e sua aceleração, gerando a curva mostrada na Figura 1(curva esta gerada pelas características do Torque aplicado).

Supondo-se que a operação, mencionada, logo acima é realizada com aceleração

constante a0, pode-se continuar a operação dentro da região denotada por A e C, porém com dificuldades na região B, desde que nesta região, o valor de as é menor do que o de a0 , desta forma, dificultando a execução desta tarefa nesta região (isto é, Ts é menor do que T0 , implicando uma diminuição no Torque).

Como conseqüência deste fato, a operação, discutida, seria seguida por as .To.

51

Figura 4: Torque Aplicado e da Aceleração angular em Problemas Não-Ideal

O objetivo, deste curso é o de estudar tais sistemas que apresentem comportamento dinâmico regular( periódico, quasi-periódico) e irregular( Caos).

A seguir, exibe-se um exemplo de um sistema vibratório, Não-ideal, com o objetivo de exibir algumas das características acima, mencionadas.

Para a análise dinâmica, deste problema, usam-se rotinas computacionais, programas em C++ e MATLAB, apresentadas, logo adiante. O modelo matemático, estudado é definido pela Figura 5 , abaixo

Figura 5- Modelo vibratório Sujeito a Excitação Paramétrica

As equações diferenciais, do movimento vibratório, definido, pela Figura 5, que são o conjunto formado pela: Equação das Vibrações Paramétricas da Haste e pela Equação que descreve o movimento da Fonte de Energia, poderão ser obtidas, fazendo-se as hipóteses:

52

• Assumiu-se que o movimento da haste tem a forma das vibrações paramétricas similares as vibrações livres.

• A mola tem uma deformação estática (fo) em adição à deformação variável:

rsinϕ-π2*y2/(4*l).

• Tomou-se:

)()(),(l

xsintytxy

∏= .

• Desprezou-se a força de inércia de um elemento da haste na direção do eixo x, a

rotação de um elemento da haste em flexão e a massa da mola;

Adotou-se a nomenclatura usual para os parâmetros físicos:

• m1 é a massa por unidade de comprimento da haste;

• Ix é o segundo momento da área de uma seção transversal da haste;

• I é o momento de inércia do rotor e;

• β é coeficiente de amortecimento linear viscoso.

Nestas condições, as equações diferenciais do movimento da haste, poder ser apresentada, como a seguir:

0)()(32 =++++ x

c ySinm

ym

y γϕωβ

&&& (*)

onde:

l

rcc

ml

clmm

EIl

PcfPP

P

m

EI

lx

x

282

1

)()1()(

12

22

41

1

21100

1

042

ππγ

ππω

−===

==−=

Note que a equação acima, descreve as vibrações paramétricas da haste, sem

porém, considerar a equação, que descreve a interação com a Fonte de Energia, no caso um motor elétrico de corrente contínua

• Equação que descreve o movimento da Fonte de Energia

53

A energia potencial V , correspondente ao modo fundamental de flexão das vibrações da haste(a energia de tensão e compressão é desprezada):

yl

Il

fcx

ESinrV

2

3

42

3

2

01)

4(

4(

2

1 ) ππϕ +−+=

A energia cinética, correpondente à ação do motor elétrico de corrente contínua é

dado por:

ϕ&2

2

1IT =

onde I é o momento de inércia do rotor. Usando as equações de Lagrange, obtém-se as equações diferenciais:

ϕϕϕϕ π Cosrl

SinrMI yfc })4

({)(2

2

01−+−= &&& (**)

onde:

• HLM −= ;

• L : é torque característico do motor e;

• H : é o momento das forças resistivas do motor Usando novas variáveis , definidas por

{ ϕϕτ && ===== 4321 ,,,, xxxxxxt }

o problema exibido na Figura 5, tem equações de movimento, na forma de estado ,

dadas por (Konokenko, 1969):

32

14333423414444

43

31241323121222

22

1

coscoscos)(,

,

xxKxsinxKxKxMKd

dxx

d

dx

xKxsinxKxKxKd

dxx

d

dx

+++==

+++==

ττ

ττ

com:

54

γ

ω

β

−=

−=

−=

−=

24

223

221

22

K

m

cK

K

Km

lI

rcK

I

rcK

I

rfcK

4

21

43

21

42

0141

π=

−=

−=

Pesquisou-se exaustivamente as soluções numéricas dos sistemas vibratório, em questão; utilizando-se vários valores dos parâmetros físicos, da fonte de energia (motor CC) e das condições iniciais, diferentes dos mencionados pelo Prof. Konokenko. O principal objetivo é analisar algumas situações topológicas não exploradas por ( Konokenko, 1969). Tomou-se os valores numéricos para os parâmetros, como sendo:

m =0.00334; c1=26; r =0.4;

Ix= 0.049; E=200000; I=0.01; l=125;

fo=0.2; ϒ=10; β=0.0001;

M=4.2- (0.0178626+g1), g1=1

com condições iniciais {x10=1;x20=0;x30=1;x40=1}, obtém-se, as Figuras 6 e 7, onde são exibidos os planos de fase x1 vs. x2 e o mapa de Poincaré associado , logo abaixo:

Figura 6: Retrato de fase x1 vs. x2

55

Figura 7: Mapa de Poincaré associado a Figura 6

O comportamento temporal da variável x1(deslocamento) é dado pela Figura 8.

Figura 8: Gráfico de x1 vs t Neste caso, usa-se o calculo da FFT, via MATLAB, para verificação do comportamento caótico, desta estrutura vibratória. Este resultado é mostrado na Figura 9

Figura 9- FFT

56

• Fenômeno do Salto(”Jump”)

O Fenômeno do Salto, foi obtido, via integração numérica, exibido pela Figura 10. Anexa-se, logo adiante, a rotina computacional que produz este resultado.

Figura 10- Resposta em freqüência( freqüência Vs. Am ) –Jump

• Rotinas Computacionais: Confecção do Mapa de Poincaré(MP)

• Rotina 2:

• A rotina em C++ que calcula os pontos para “plotar” o mapa de Poincaré, foi criada baseada na seguinte lógica, descrita no algoritmo, descrito por(Parker e Chuá, 1998).

• programa para calcular o (MP), primeiro integra as equações do sistema

vibratório, através do integrador RK4, e em cada dupla de pontos provenientes desta integração é chamada a rotina Poincaré onde se verifica se houve” furo” no plano escolhido que irá mapear a trajetória no espaço.

• Essa verificação é feita no passo 3 do algoritmo e em caso positivo é feito a

interpolação e, em seguida verifica-se se o valor está dentro de uma faixa de controle para visualização (artifício para visualizar pontos tridimensionais em duas dimensões.

• Tendo estes pontos, eles são armazenados em um arquivo em Formato ASCII e,

o programa prossegue processando os outros pontos no intervalo de integração especificado.

57

• Terminado o programa, é executada uma rotina, no MATLAB, que lê esses pontos do arquivo, disponibilizando-os na memória podendo, ser então, “ plotado” o mapa de Poincaré(MP) quando se deseja, através do comando plot do MATLAB.

• A rotina para o MP escrita, nesta linguagem, tem a mesma lógica, mas não

tem um integrador imbutido nela, o que torna necessário, primeiro efetuar a integração numérica, usando um integrador numérico, no caso RK4 que foi escrito em linguagem de programação C++ armazenando os ponto em arquivo no formato ASCII .

Os Programas usados, na análise do problema, em tela, foram elaborados, a partir

dos algoritmos , logo acima, discutidos e apresentados, a seguir.

Programas feito no MATLAB %**********************************************************************% %Programa 1: Carregar % % Programa para carregar um arquivo ascii (de pontos), que é solução de um sistema de equa- % %ções diferenciais, e plotar o gráfico destes mesmos %%********************************************************************% load (caminho)nomearq.dat t=nomearq(:,1); x1=nomearq(:,2); %CARREGANDO EM VARIÁVEIS AS COLUNAS DO ARQUIVO x2=nomearq(:,3); x3=nomearq(:,4); x4=nomearq(:,5); %Plotando o gráfico t vs x1 plot(t,x1) pause % Plotando o gráfico t vs x2

plot(t,x2) pause % Plotando o gráfico t vs x3 plot(t,x3) pause % Plotando o gráfico t vs x4 plot(t,x4) pause % Plotando o plano de fase x1 vs x2 plot(x1,x2)

//**************************************************************// // Programa 2: Integrador Numérico // //Comentário: // // Programa para integrar um sistema de equações diferenciais, // // devendo este ser um sistema de primeira ordem. As saídas (solução): // // tempo,x1,x2,x3,x4 são salvas em arquivo ascii SOL.DAT para serem usadas // // posteriormente // // Este programa usa integrador de Quarta ordem da família Runge-Kutta // //**************************************************************// #include<stdio.h>#include<conio.h>#include<math.h>#include <stdlib.h>#include <process.h>

//Definição da constante Pi #define PI 3.1415

58

//Declaração das variáveisdouble a,b,n,h,t;double ci[5], w[5], k1[5], k2[5], k3[5], k4[5],acumula[5],xp[5],normal[5];unsigned long i,j,equa,x,saida;FILE *fptr,*ind; //Declaração das constantesdouble m =0.0034; double c1 =26; double r =0.4; double Ix =0.049; double E =200000; double I =0.01; double l =125; double g1 =1; double fo =0.2; double beta=0.0001; double P1 =powl(PI/l,2)*E*Ix; double Po =fo*c1; double m1=2*m/l; double c2 =-((powl(PI,2)*r*c1) /(2*l)) ; double gama=powl(PI,4)*c1/(8*l*l*m); double omega2=powl(PI/l,4)*E*(Ix/m1)*(1-(Po/P1)); double k21 = -omega2; double k22 = -beta/m; double k23 = -c2/m; double k24 = -gama; double k41 = -(c1*fo*r)/I; double k42= -(c1*r*r)/I; double k43= c1*(powl(PI,2)*r)/(4*l*I); double k44= 1/I; // Definição das equações diferenciais em 1º ordem do problema estudado.double f(int indice,double t,double x1,double x2,double x3,double x4){ double resp,aux; //sistema de equacoes diferencias de 1 ordem if(indice==1) resp = x2; if(indice==2) resp = k22*x2 + k21*x1 + k23*sin(x3)*x1 +k24*powl(x1,3); if(indice==3) resp = x4; if(indice==4) resp = k44*(4.2-(0.0178626+g1)*x4) + +k41*cos(x3)+k42*sin(x3)*cos(x3)+ +k43*powl(x1,2)*cos(x3); return(resp);

} void main(void) { equa =4; a =0; b =10; n =3000; h = (b-a)/n; t = a; // condicoes iniciais de x1,x2,x3,x4, respectivamente. ci[1] =1; ci[2] =0; ci[3] =1; ci[4] =1; for(j=1;j<=equa;j++) w[j] = ci[j]; printf("\n%f %f",t,w[1]); //Abrindo o arquivo que conterá a solução fptr = fopen("C:\\Osmar\\pemc\\sol.dat","w+"); fprintf(fptr,"%f %f %f %f %f\n",t,w[1],w[2],w[3],w[4]); for(i=1;i<=n;i++) { for(j=1;j<=equa;j++) k1[j] = h * f( j, t, w[1], w[2], w[3], w[4]); for(j=1;j<=equa;j++) k2[j] = h * f( j, t+h/2, w[1]+k1[1]/2, w[2]+k1[2]/2, w[3]+k1[3]/2, w[4]+k1[4]/2); for(j=1;j<=equa;j++) k3[j] = h * f( j, t+h/2, w[1]+k2[1]/2, w[2]+k2[2]/2, w[3]+k2[3]/2, w[4]+k2[4]/2); for(j=1;j<=equa;j++) k4[j] = h * f( j, t+h, w[1]+k3[1], w[2]+k3[2], w[3]+k3[3], w[4]+k3[4]); for(j=1;j<=equa;j++) w[j] = w[j] + (k1[j] + 2*k2[j] + 2*k3[j] + k4[j])/6; t = a + i*h; //comando para ver em tela a solução encontrada printf("\n%f %f ",t,w[1]); fprintf(fptr,"%f %f %f %f %f\n",t,w[1],w[2],w[3],w[4]); } //Fechando o arquivo fclose(fptr);getch();}

59

//**************************************************************// // Programa 3: Mapa de Poincaré // //Programa em C, que tem como saída, pontos em um arquivo ascii // // POINCAR.DAT para elaboração do Mapa de Poincaré. O integrador (RK4) // //fornece quatro saídas e a função poincare faz a projeção destas saídas em // //em plano tridimensional, podendo posteriormente ser plotado o mapa. // //**************************************************************// #include<stdio.h> #include<conio.h> #include<math.h> #include <stdlib.h> #include <process.h> //Definição da constante π #define PI 3.1415 //Declaração das variáveis double a,b,n,h,t; double ci[5], w[5], k1[5], k2[5], k3[5], k4[5],acumula[5],xp[5],normal[5]; unsigned long i,j,equa,x,saida; FILE *fptr,*ind; // Definicao das constantes double m =0.00334; double c1 = 26; double r = 0.2; double Ix =0.049; double E =200000; double I =0.01; double l =125; double g1 =1; double fo =0.2; double beta=0.0001; double P1 =powl(PI/l,2)*E*Ix; double Po =fo*c1; double c2 =-((powl(PI,2)*r*c1) /(2*l)) ; double gama=powl(PI,4)*c1/(8*l*l*m); double omega2=powl(PI/l,4)*E*(Ix/m)*(1-(Po/P1)); double k21 = -omega2; double k22 = -beta/m; double k23 = -c2/m; double k24 = -gama; double k41 = -(c1*fo*r)/I; double k42= -(c1*r*r)/I; double k43= c1*(powl(PI,2)*r)/(4*l*I); double k44= 1/I; //rotina para calcular o Mapa de poincare

void poincare() { int o,k; double alf1=0; double alf2=0; double xmed1[5]; double xmed2[5]; double xmed[5]; for (o=1;o<=equa;o++) { alf1=alf1+(normal[o] * (acumula[o]-xp[o])); } for (o=1;o<=equa;o++) { alf2=alf2+normal[o]*(w[o]-xp[o]); } if (alf1>0) { if (alf1*alf2<0) { saida=saida+1; for(o=1;o<=equa;o++) { xmed1[o]=(alf2/(alf2-alf1))*acumula[o]; xmed2[o]=(alf1/(alf1-alf2))*w[o]; xmed[o]=xmed1[o]+xmed2[o]; } fprintf(ind,"%f %f %f %f\n",xmed[1],xmed[2],xmed[3],xmed[4]); } } if(saida>200000) exit(0); } // Definicao das equacoes diferenciais do problema estudado. double f(int indice,double t,double x1,double x2,double x3,double x4) {

60

double resp; //Sistema de equações diferencias de 1º ordem if(indice==1) resp = x2; if(indice==2) resp = k22*x2 + k21*x1 + k23*sin(x3)*x1 +k24*powl(x1,3); if(indice==3) resp = x4; if(indice==4) resp = k44*(4.2-(0.0178626+g1)*x4) + +k41*cos(x3)+k42*sin(x3)*cos(x3)+ +k43*powl(x1,2)*cos(x3); return(resp); } void main(void) { clrscr(); ldiv_t x; //Declaração das constantes equa =4; a =0; b =10000000; n =1000000000; h = (b-a)/n; t = a; // Condições iniciais de x1,x2,x3,x4, respectivamente. ci[1] =1; ci[2] =0; ci[3] =1; ci[4] =1; //Tomando um ponto no plano desejado (xp) e declarando seu vetor normal // ao plano escolhido xp[1]=-0.2; xp[2]=0; xp[3]=0; xp[4]=3.5; normal[1]=0; normal[2]=0; normal[3]=0; normal[4]=1; for(j=1;j<=equa;j++) w[j] = ci[j]; //Abrindo o arquivo para escrita ind = fopen("C:\\Osmar\\pemc\\poincar.dat","w+"); for(i=1;i<=n;i++) { for (j=1;j<=equa;j++) acumula[j]=w[j];

for(j=1;j<=equa;j++) k1[j] = h * f( j, t, w[1], w[2], w[3], w[4]); for(j=1;j<=equa;j++) k2[j] = h * f( j, t+h/2, w[1]+k1[1]/2, w[2]+k1[2]/2, w[3]+k1[3]/2, w[4]+k1[4]/2); for(j=1;j<=equa;j++) k3[j] = h * f( j, t+h/2, w[1]+k2[1]/2, w[2]+k2[2]/2, w[3]+k2[3]/2, w[4]+k2[4]/2); for(j=1;j<=equa;j++) k4[j] = h * f( j, t+h, w[1]+k3[1], w[2]+k3[2], w[3]+k3[3], w[4]+k3[4]); for(j=1;j<=equa;j++) w[j] = w[j] + (k1[j] + 2*k2[j] + 2*k3[j] + k4[j])/6; t = a + i*h; x=ldiv(i,2); if (x.rem==0) poincare(); } //Fechando o arquivo fclose(ind); getch(); }

61

//**************************************************************// // Programa 4 : Resposta em Frêquencia // //Programa em C que obtém, através do integrador RK4, valores da amplitude // //e da frequência em um determinado intervalo, variando E1, e armazena-os // // em um arquivo de nome RESP.DAT. Com este resultado, torna-se possível a// //construção um gráfico para visualizar a resposta em frequência. As saídas // // armazenadas são ml(curva característica do motor),amplitude e frequência //// salvas em arquivo ascii para serem usadas posteriormente //// //**************************************************************// #include<stdio.h> #include<conio.h> #include<math.h> //Definição da constante π #define PI 3.1415 //Definição das variáveisdouble a,b,ml1,ml,n,h,t; double ci[5], w[5], k1[5], k2[5], k3[5], k4[5]; long int i,j,k,equa,count; FILE *fptr; //Definição das constantes double m =0.0034; double c1 =26; double r = 0.4;double Ix =0.049; double E =200000; double I =0.01; double l =125; double g1 =1; double fo =0.2; double beta=0.001; double P1 =powl(PI/l,2)*E*Ix; double Po =fo*c1; double m1=2*m/l; double c2 =-((powl(PI,2)*r*c1) /(2*l)) ; double gama=powl(PI,4)*c1/(8*l*l*m); double omega2=powl(PI/l,4)*E*(Ix/m1)*(1-(Po/P1)); double k21 = -omega2; double k22 = -beta/m; double k23 = -c2/m; double k24 = -gama; double k41 = -(c1*fo*r)/I;

double k42= -(c1*r*r)/I; double k43= c1*(powl(PI,2)*r)/(4*l*I); double k44= 1/I; // Definição das equações diferenciais do problema estudado. double f(int indice,double t,double x1,double x2,double x3,double x4) { double resp; //double ml=ml1*(1-0.00425233*x4); ml=ml1*(1 -0.00425233*x4); //sistema de equações diferencias de 1 ordem if(indice==1) resp = x2; if(indice==2) resp = k22*x2 + k21*x1 + k23*sin(x3)*x1 +k24*powl(x1,3); if(indice==3) resp = x4; if(indice==4) resp = k44*(ml-(g1*x4)) + +k41*cos(x3)+k42*sin(x3)*cos(x3)+ +k43*powl(x1,2)*cos(x3); return(resp); } void main(void) { //Declaração de variáveis double tempo,tempo_max,vw,maximo; double t1,t2,t3,t4; double df1,df2,df3; double periodo,frequencia,amplitude; clrscr();

62

// Definição das constantes equa=4; a = 0; b =30; n =11400; h = (b-a)/n; t = a; fptr = fopen("C:\\osmar\\pemc\\resp.dat","w+"); for(k=0;k<=40;k++) { ml1 =2 + k*0.021; // condições iniciais de x1,x2,x3,x4, respectivamente. ci[1] =1; ci[2] =0; ci[3] =1; ci[4] = 1; for(j=1;j<=equa;j++) w[j] = ci[j]; tempo = 0; tempo_max = 0; maximo = 0; t = 0; t1=t2=t3=t4=0; for(i=1;i<=n;i++) { for(j=1;j<=equa;j++) k1[j] = h * f( j, t, w[1], w[2], w[3], w[4]); for(j=1;j<=equa;j++) k2[j] = h * f( j, t+h/2, w[1]+k1[1]/2, w[2]+k1[2]/2,w[3]+k1[3]/2, w[4]+k1[4]/2); for(j=1;j<=equa;j++) k3[j] = h * f( j, t+h/2, w[1]+k2[1]/2, w[2]+k2[2]/2, w[3]+k2[3]/2, w[4]+k2[4]/2); for(j=1;j<=equa;j++) k4[j] = h * f( j, t+h, w[1]+k3[1], w[2]+k3[2],

w[3]+k3[3], w[4]+k3[4]); for(j=1;j<=equa;j++) w[j] = w[j] + (k1[j] + 2*k2[j] + 2*k3[j] + k4[j])/6; t = a + i*h; if(w[1]>=0) { //Encontrando a amplitude if((vw < maximo) && (maximo > w[1])) { t4=t3; t3=t2; t2=t1; t1 = tempo_max; amplitude = maximo; } } vw = maximo; maximo = w[1]; tempo = tempo_max; tempo_max = t; }//fim do for ate n //Calcula a frequencia df1 = t1-t2; df2 = t2-t3; df3 = t3-t4; periodo = (df1+df2+df3)/3; frequencia = 2*acos(-1)/periodo; //saida para o arquivo fprintf(fptr,"%f %f %f\n",ml1, frequencia,amplitude); //saida para a tela printf("\n %f %f %f",ml1,amplitude,frequencia); } fclose(fptr); getch(); }

63

Finalmente, toma-se o modelo, de problema não-ideal, ilustrado pela Figura 11

Figura 11- Problema não-ideal

O problema dado pela Figura 11, consiste: • De uma massa m1, de uma mola de coeficiente de rigidez( linear) k2 e

coeficiente de amortecimento c1 . • No corpo de m1, coloca-se uma fonte não-ideal de potência limitada( Motor CC)

• O motor CC tem rotor de momento de inércia J e massa excêntrica m0

situada à uma r do eixo de rotação do .

• Através de uma mola linear, com rigidez k2 e coeficiente de amortecimento c2, toma-se um corpo de massa m2 atachado à massa m1.

As equações diferenciais do movimento são:

( ) ( )m x k x c x k x x c x x m r m r1 1 1 1 1 1 2 2 1 2 2 1 0

2

0&& & & & cos & sin ,= − − + − + − + +ω ϕ ω ϕ

( ) ( )m x k x x c x x2 2 2 2 1 2 2 1&& & &= − − − −

( )J L H m x& && sinω ω ϕ= − + 0 1 ( )&L aL b U= − − +ω ω ,

onde:

• L é o torque gerado pelo motor, • ( )H ω é o torque resistivo do motor,

• a e b são constantes que dependem do tipo de motor elétrico de corrente contínua tomado para estudo e,

• ( )U ω é a tensão do motor cc.

64

O objetivo deste estudo é o de sintetizar um controle ótimo de Passagem pela ressonância, através do primeiro pico, na região de anti-ressonância, usando o método de regularização de Tikhonov’s(Tikhnov e Arsenin, 1977).

Usando variáveis adimensionais, e usando o método de Tikhonov’s obteve-se os resultados, ilustrados, através das Figuras 12 e 13 logo abaixo(Cheshankov et al., 1999):

Time

0 100 200 300 400 500 600 700

Dis

pla

cem

ent

-2

-1

0

1

2

Figure 1 2 : Passagem pela Ressonância, sem Controle

Time

0 100 200 300 400 500 600 700

Dis

pla

cem

ent

-2

-1

0

1

2

Figure 13: Passagem pela Ressonância, sem Controle

65

66

• Cheshankov, B. I. , D.T. Ruschev, J.M. Balthazar, H.I. Weber. A Synthesis of an Optimal Control of a Nonideal Energy Source by Use of Tichonov’s Method of Regularization. IN: ICONNE 2000, to be held in Campos do Jordão, July 31- August4, 2000, 2000, aceito para publicação.