19
Departamento de Engenharia Mecânica Desenvolvimento de uma Metodologia para a Implementação de Curvas Cônicas na Elaboração de Trajetórias Ótimas de um Veículo Aluno: Ramon Felipe Brandão do Nascimento Orientador: Mauro Speranza Neto Introdução Uma das maiores dificuldades da engenharia de corridas é a descoberta do melhor traçado em uma pista para que um veículo possa percorrê-la em menor tempo. Um dos fatos que tornam essa tarefa difícil é a concatenação de traçados ótimos isolados em diferentes curvas consecutivas, uma vez que na maioria das vezes as trajetórias estipuladas nas diferentes curvas não são contínuas. Isso se deve aos diferentes formatos de trajetórias possíveis, visto que a trajetória ótima para uma curva deve estar diretamente relacionada às curvas seguintes. Diante dessa dificuldade, para que seja possível obter o melhor percurso de um veículo em um circuito, podem ser feitas algumas aproximações quanto as curvas que irão unir os melhores traçados isolados. Isso se deve ao fato de, muitas vezes, essas curvas possuírem modelos difíceis de serem trabalhados, como a Clotóide, por exemplo, que é um é uma curva de curvatura variável. A partir da tentativa de se obter esses traçados, percebeu-se que para facilitar esse processo, algumas curvas complexas poderiam ser substituídas por curvas cônicas (parábolas, hipérboles e elipses), tentando obter assim uma boa aproximação. O motivo de serem escolhidas cônicas é o fato delas possuírem seus parâmetros bem definidos, o que permite, após a elaboração de um modelo, obter analiticamente todas as variáveis do veículo na curva, como suas acelerações normal e tangencial, velocidades e acelerações em relação aos referenciais globais. Objetivo Uma vez que a tarefa de se encontrar trajetórias ótimas não é um processo muito simples, esse trabalho se concentrou no estudo das curvas cônicas e na elaboração de um programa em “Matlab” que reconhece a curva em sua forma paramétrica e fornece todos os parâmetros de uma partícula orientada que a utiliza como trajetória. Entretanto, tem-se como objetivo a continuação desse trabalho no próximo ano para que, com isso, possam ser feitas análises de traçados ideais de pistas reais a partir de trechos de curvas cônicas.

Desenvolvimento de uma Metodologia para a Implementação de ... · Desenvolvimento de uma Metodologia para a Implementação de Curvas Cônicas na Elaboração de Trajetórias Ótimas

  • Upload
    others

  • View
    7

  • Download
    0

Embed Size (px)

Citation preview

Departamento de Engenharia Mecânica

Desenvolvimento de uma Metodologia para a Implementação de Curvas

Cônicas na Elaboração de Trajetórias Ótimas de um Veículo

Aluno: Ramon Felipe Brandão do Nascimento

Orientador: Mauro Speranza Neto

Introdução

Uma das maiores dificuldades da engenharia de corridas é a descoberta do melhor

traçado em uma pista para que um veículo possa percorrê-la em menor tempo. Um dos fatos

que tornam essa tarefa difícil é a concatenação de traçados ótimos isolados em diferentes curvas

consecutivas, uma vez que na maioria das vezes as trajetórias estipuladas nas diferentes curvas

não são contínuas. Isso se deve aos diferentes formatos de trajetórias possíveis, visto que a

trajetória ótima para uma curva deve estar diretamente relacionada às curvas seguintes.

Diante dessa dificuldade, para que seja possível obter o melhor percurso de um veículo

em um circuito, podem ser feitas algumas aproximações quanto as curvas que irão unir os

melhores traçados isolados. Isso se deve ao fato de, muitas vezes, essas curvas possuírem

modelos difíceis de serem trabalhados, como a Clotóide, por exemplo, que é um é uma curva

de curvatura variável.

A partir da tentativa de se obter esses traçados, percebeu-se que para facilitar esse

processo, algumas curvas complexas poderiam ser substituídas por curvas cônicas (parábolas,

hipérboles e elipses), tentando obter assim uma boa aproximação. O motivo de serem

escolhidas cônicas é o fato delas possuírem seus parâmetros bem definidos, o que permite, após

a elaboração de um modelo, obter analiticamente todas as variáveis do veículo na curva, como

suas acelerações normal e tangencial, velocidades e acelerações em relação aos referenciais

globais.

Objetivo

Uma vez que a tarefa de se encontrar trajetórias ótimas não é um processo muito

simples, esse trabalho se concentrou no estudo das curvas cônicas e na elaboração de um

programa em “Matlab” que reconhece a curva em sua forma paramétrica e fornece todos os

parâmetros de uma partícula orientada que a utiliza como trajetória. Entretanto, tem-se como

objetivo a continuação desse trabalho no próximo ano para que, com isso, possam ser feitas

análises de traçados ideais de pistas reais a partir de trechos de curvas cônicas.

Departamento de Engenharia Mecânica

Metodologia

Foi desenvolvido um modelo para cada curva cônica separadamente, os quais serão

analisados a seguir:

1) Círculo

Dentre as curvas analisadas, o círculo é o modelo mais simples, uma vez que seu raio de

curvatura permanece constante ao longo da trajetória. Por esse fato, iniciar-se-á a análise por

ele.

A equação cartesiana característica de uma curva cônica é dada pela expressão:

A partir dessa equação, é possível descobrir qual das cônicas ela descreve a partir das

relações:

{𝐵2 = 4𝐴𝐶 → 𝑃𝑎𝑟á𝑏𝑜𝑙𝑎𝐵2 > 4𝐴𝐶 → 𝐻𝑖𝑝é𝑟𝑏𝑜𝑙𝑒

𝐵2 < 4𝐴𝐶 → 𝐶í𝑟𝑐𝑢𝑙𝑜 𝑜𝑢 𝐸𝑙𝑖𝑝𝑠𝑒

Para o círculo, além de 𝐵2 < 4𝐴𝐶, o valor de B deve ser nulo e A=C.

A partir de um sistema de mudança de coordenadas de translação e exercícios de álgebra é

possível se determinar os parâmetros característicos do círculo. A posição de seu centro no eixo

X e Y é dada por 𝑋0 =−𝐷

2𝐴 e 𝑌0 =

−𝐸

2𝐴. Além disso, o valor do Raio de Curvatura é dado pela

expressão √−𝐹′

𝐴, onde 𝐹′ é dado por 𝐴(𝑋0 + 𝑌0) + 𝐷𝑋0 + 𝐸𝑌0 + 𝐹. Com isso a mudança de

coordenadas a equação cartesiana fica 𝑋′2 + 𝑌′2 = |𝐹′

𝐴|.

Simulação

Para se analisar o círculo como a trajetória de uma partícula orientada são dados como

parâmetros iniciais sua Velocidade Inicial (𝑉0), a variação de sua Aceleração Tangencial ao

longo do tempo, permanecendo constante a cada intervalo de tempo até que assuma outro valor,

e a sua respectiva equação cartesiana.

A parametrização do círculo utilizada é dada por 𝛼(𝜃) = (𝑅𝑐𝑜𝑠(𝜃), 𝑅𝑠𝑖𝑛(𝜃)), fornecendo

assim as coordenadas X e Y do círculo nas coordenadas locais. Utilizando a equação de

Movimento Circular Uniformemente Acelerado 𝜃(𝑡) = 𝜃0 + 𝜔0𝑡 +𝛼𝑡2

2, no qual 𝛼 =

𝐴𝑡

𝑅 e

𝜔0 = 𝑉0

𝑅, é possível determinar a posição da partícula em função do tempo. Além disso, ao

derivar a equação 𝑉.𝛼(𝜃)

𝑅, encontra-se 𝑉(𝑡) = (−𝑉𝑠𝑖𝑛(𝑡), 𝑉𝑐𝑜𝑠(𝑡)), onde 𝑉 = 𝑉0 + 𝐴𝑡𝑡.

Depois, ao aplicar o sistema de transformação de coordenadas de translação {𝑋 = 𝑋′ + 𝑋0

𝑌 = 𝑌′ + 𝑌0 ,

determina-se a posição do círculo no sistema de coordenadas global.

Departamento de Engenharia Mecânica

Abaixo segue um exemplo de uma simulação em “Matlab” de um círculo utilizando a

equação cartesiana 10𝑋2 + 10𝑌2 + 5𝑋 + 7𝑌 + 100.

clear; clc;

%Coeficientes

A=10;

B=0;

C=10;

D=5;

E=7;

F=100;

%Variaveis Iniciais

V0=0; %Velocidade Inicial

(m/s)

X0= -D/2*A; %Posição do

centro no eixo X

Y0= -E/2*A; %Posição do

centro no eixo Y

Fl = A*X0^2 + A*Y0^2 + D*X0

+ E*Y0 + F;

R= sqrt(Fl/A); %Raio de

Curvatura

ti=0; Tempo Inicial (s)

theta0=0; %Angulo Inicial

(rad)

w0=V0/R; %Velocidade Angular

Inicial

theta(1)=theta0;

w(1)=w0;

At(1)=1;

%Aceleracao

Tangencial

(m/s^2)

Xl(1)= R;

Yl(1)=0;

Vx(1)=0;

Vy(1)=V0;

V(1)=V0;

X(1)= X0 + Xl(1);

Y(1) = Y0 + Yl(1);

An(1)= V(1)^2/R; %Aceleracao

Normal (m/s^2)

t(1)=ti; %Tempo Inicial

tf = 50; %Tempo Final (s)

dt = .01; %Intervalo de

Simulacao

for i=1:tf/dt

t(i+1)= i*dt;

if t(i+1)<20

At(i+1) = 1;

else

At(i+1) = 0;

end

V(i+1)=V(i) + At(i)*dt;

alfa = At(i)/R;

w(i+1)=V(i+1)/R;

if At(i)==0

theta(i+1)=theta(i)

+ w(i)*dt;

theta0=theta(i+1);

w0=w(i+1);

ti = t(i+1);

else

theta(i+1)=theta0 +

w0*(t(i+1)-ti) +

alfa/2*(t(i+1)-ti)^2;

end

An(i+1) = (V(i+1)^2)/R;

%m/s^2

Xl(i+1) =

R*cos(theta(i+1));

Yl(i+1) =

R*sin(theta(i+1));

Vx(i+1) = -

V(i+1)*sin(theta(i+1));

Vy(i+1)=V(i+1)*cos(theta(i+1

));

X(i+1)= Xl(i+1) + X0;

Y(i+1) = Yl(i+1)+ Y0;

End

Departamento de Engenharia Mecânica

2) Elipse

As próximas análises serão mais difíceis pelo fato do Raio de Curvatura variar em cada

ponto da trajetória, diferentemente do círculo. Para que ele pudesse ser descoberto mais

facilmente, foi desenvolvido um programa em “Maple” que, ao ser fornecida a equação

paramétrica de qualquer curva, é dado o módulo do Raio de Curvatura e o Vetor Normal,

facilitando assim a descoberta do centro instantâneo. Segue abaixo o programa exemplificado

para a Elipse:

>

Derivada de r pela variável 𝜃

Derivada do arco S da curva pela variável 𝜃

Departamento de Engenharia Mecânica

`

Vetor Tangente

Segunda derivada de r pela variável 𝜃

Segunda derivada do arco S da curva pela variável 𝜃

Raio de Curvatura

Vetor Normal

Departamento de Engenharia Mecânica

Para se descobrir o centro de curvatura instantâneo, basta fazer �̅� + 𝜌𝑒𝑛̅̅ ̅, onde �̅� é a curva

parametrizada, 𝜌 é o raio de curvatura e 𝑒𝑛̅̅ ̅ o vetor normal a curva

Utilizando sistemas de transformação de coordenadas de translação e rotação,

respectivamente, e utilizando alguns artifícios algébricos na equação cartesiana da elipse é

possível determinar seus parâmetros e sua equação em novas coordenadas. A posição de seu

centro nas coordenadas X e Y é dada por e .

Além disso, o ângulo de rotação utilizado na transformação é dado por

ou , caso A=C. Sua equação em coordenadas locais fica 𝑋′′2

𝑎+

𝑌′′2

𝑏= 1, onde 𝑎 =

𝐹′

𝑢 e 𝑏 =

𝐹′

𝑣. 𝐹′ = 𝐴𝑋0 + 𝐵𝑋0𝑌0 + 𝐶𝑌0 + 𝐷𝑋0 + 𝐸𝑌0 + 𝐹, 𝑢 =

𝐴𝑐𝑜𝑠2(𝜑) + 𝐵𝑠𝑖𝑛(𝜑) cos(𝜑) + 𝐶𝑠𝑖𝑛2(𝜑) e 𝑣 = 𝐴𝑠𝑖𝑛2(𝜑) − 𝐵𝑠𝑖𝑛(𝜑) cos(𝜑) + 𝐶𝑐𝑜𝑠2(𝜑).

Simulação

Os parâmetros iniciais de simulação são os mesmos fornecidos no caso do círculo.

Entretanto, como o Raio de Curvatura varia ao longo da trajetória, o algoritmo utilizado nesse

caso será diferente do caso do círculo. Todavia, há o benefício desse modelo poder ser

reutilizado no modelo das outras curvas analisadas. O diagrama a seguir representa o modelo

que foi implementado em “Matlab”:

Onde:

𝐴𝑡 → Aceleração Tangencial

Departamento de Engenharia Mecânica

𝛼 → Aceleração Rotacional

𝐴𝑛 → Aceleração Normal

𝑉 𝑒 𝑉0 → Velocidades Instantânea e Inicial

𝜌 → Raio de Curvatura Instantâneo

𝜔 → Velocidade Angular

𝜃 → Posição Angular

𝑋′′𝑒 𝑌′′ → Coordenadas no Referencial Local

𝑉′′𝑥 𝑒 𝑉′′𝑦 → Velocidades no Referencial Local

Abaixo segue um exemplo de uma simulação em “Matlab” de uma elipse utilizando a

equação cartesiana 90𝑋2 + 30𝑋𝑌 + 60𝑌2 − 4000𝑋 − 3000𝑌 + 7000:

clear; clc;

%Coeficientes

A=90;

B=30;

C=60;

D=-4000;

E=-3000;

F=7000;

%Variaveis Iniciais

V0=3; %Velocidade Inicial em m/s

V(1)=V0;

X0= -((-B * E + 2 * C * D + B * D)

/ A / (2 * C + B)) /

0.2e1;%Posição do centro no eixo X

Y0= -E / (2 * C + B);%Posição do

centro no eixo Y

Fl = A*X0^2 + B*X0*Y0 + C*Y0^2 +

D*X0 + E*Y0 + F;

if A==C

psi=pi/4;

else

psi= atan((B/(A-

C)))*(1/2);%Angulo de Rotacao

end

u = A*(cos(psi))^2 +

B*sin(psi)*cos(psi) +

C*(sin(psi))^2;

v = A*(sin(psi))^2 -

B*sin(psi)*cos(psi) +

C*(cos(psi))^2;

a= sqrt(-Fl/u);%Coeficiente de Xll

b= sqrt(-Fl/v);%Coeficiente de Yll

theta0=0;%Angulo Inicial

theta(1)=theta0;

w0=0;%Velocidade Angular Inicial

w(1)=w0;

At(1)=1;%Aceleracao Tangencial

Inicial

%Posicoes em direfentes

coordenadas

Xll(1)= a;

Yll(1)=0;

Xl(1) = cos(psi)*Xll(1) -

sin(psi)*Yll(1);

Yl(1) = sin(psi)*Xll(1) +

cos(psi)*Yll(1);

X(1)= X0 + Xl(1);

Y(1) = Y0 + Yl(1);

R(1) = RdCE(0, a, b);%Raio de

Curvatura Instantaneo

N= VNE(0,a,b);%Vetor Normal

NXll(1) = R(1)*N(1,1)+a;

NYll(1) = R(1)*N(2,1);

NXl(1) = cos(psi)*NXll(1) -

sin(psi)*NYll(1);

NYl(1) = sin(psi)*NXll(1) +

cos(psi)*NYll(1);

NX(1)= X0 + NXl(1);

NY(1) = Y0 + NYl(1);

%Velocidades em Diferentes

Coordenadas

Vxll(1) = -

V0*a*sin(theta(1))/(sqrt((a*sin(th

eta(1)))^2 +

(b*cos(theta(1)))^2));

Vyll(1) =

V0*b*cos(theta(1))/(sqrt((a*sin(th

eta(1)))^2 +

(b*cos(theta(1)))^2));

Vx(1)=cos(psi)*Vxll(1) -

sin(psi)*Vyll(1);

Departamento de Engenharia Mecânica

Vy(1)=sin(psi)*Vxll(1) +

cos(psi)*Vyll(1);

An(1)= 0;%Aceleracao Normal

Inicial

t(1)=0;%Tempo Inicial

A(1)=0;%Modulo da Aceleracao

Inicial

tf = 25; %Tempo Final

dt = .005; %Intervalo de Simulacao

for i=1:1:tf/dt

t(i+1)=i*dt;

if t(i+1)<5

At(i+1) = 1;

elseif t(i+1)>=5 && t(i+1)<7

At(i+1) = .5;

else

At(i+1) = 0;

end

V(i+1)=V(i)+At(i)*dt;

An(i+1)=(V(i))^2/R(i);

A(i+1)=sqrt((An(i+1))^2 +

(At(i+1)^2));

w(i+1)=V(i)/R(i);

theta(i+1)=theta(i)+w(i)*dt +

At(i+1)/R(i)*dt^2;

R(i+1) = RdCE(theta(i+1), a,

b);

Xll(i+1) = a*cos(theta(i+1));

Yll(i+1) = b*sin(theta(i+1));

Xl(i+1) = cos(psi)*Xll(i+1) -

sin(psi)*Yll(i+1);

Yl(i+1) = sin(psi)*Xll(i+1) +

cos(psi)*Yll(i+1);

X(i+1)= X0 + Xl(i+1);

Y(i+1) = Y0 + Yl(i+1);

N = VNE(theta(i+1), a, b);

NXll(i+1) = R(i+1)*N(1,1)+

Xll(i+1);

NYll(i+1) = R(i+1)*N(2,1)+

Yll(i+1);

NXl(i+1) = cos(psi)*NXll(i+1)

- sin(psi)*NYll(i+1);

NYl(i+1) = sin(psi)*NXll(i+1)

+ cos(psi)*NYll(i+1);

NX(i+1)= X0 + NXl(i+1);

NY(i+1) = Y0 + NYl(i+1);

Vxll(i+1) = -

V(i+1)*a*sin(theta(i+1))/(sqrt((a*

sin(theta(i+1)))^2 +

(b*cos(theta(i+1)))^2));

Vyll(i+1)=V(i+1)*b*cos(theta(i+1))

/(sqrt((a*sin(theta(i+1)))^2 +

(b*cos(theta(i+1)))^2));

Vx(i+1)=cos(psi)*Vxll(i+1) -

sin(psi)*Vyll(i+1);

Vy(i+1)=sin(psi)*Vxll(i+1) +

cos(psi)*Vyll(i+1);

end

Departamento de Engenharia Mecânica

3) Parábola

A Parábola, na sua forma simples com seus eixos característicos locais, pode possuir

simetria no eixo 𝑋′′ ou 𝑌′′ e sua equação assume a forma 𝑌′′2 = 2𝑝𝑋′′ ou 𝑋′′2 = 2𝑝𝑌′′, respectivamente, sendo p o dobro da distância focal da parábola. Se p for positivo, a parábola

pode possuir concavidade para cima ou para a direita e, se for negativo, para baixo ou para a

esquerda. Por esse fato, a parametrização da parábola por coordenadas polares pode assumir

quatro formas distintas, dependendo do valor de p e do seu eixo de simetria. Após serem

introduzidos outros parâmetros dessa curva, serão apresentadas as quatro parametrizações

possíveis.

Aplicando-se o sistema de transformações de rotação e de translação, nessa ordem, na

equação cartesiana da parábola e desenvolvendo alguns cálculos é possível obter a posição do

vértice da parábola (𝑋0 𝑒 𝑌0), o valor de p e descobrir dentre u e v qual é zero, determinando

assim qual o eixo de simetria da parábola.

Caso 𝑢 ≅ 0 o eixo de simetria é 𝑋′′, 𝑝 =−𝑘

2𝑣, 𝑋0 =

−𝐹

𝑘+

𝑗2

4𝑣𝑘 e 𝑌0 =

−𝑗

2𝑣.

1. Se 𝑝 > 0, sua equação polar é 𝜌 =𝑝

1−cos (𝜃).

2. Se 𝑝 < 0, sua equação polar é 𝜌 =𝑝

1+cos (𝜃).

Caso 𝑣 ≅ 0 o eixo de simetria é 𝑌′′, 𝑝 =−𝑗

2𝑢, 𝑋0 =

−𝑘

2𝑢 e 𝑌0 =

−𝐹

𝑗+

𝑘2

4𝑢𝑗

3. Se 𝑝 > 0, sua equação polar é 𝜌 =𝑝

1−sin(𝜃).

4. Se 𝑝 < 0, sua equação polar é 𝜌 =𝑝

1+sin(𝜃).

U e v assumem o mesmo valor que na elipse, 𝑘 = 𝐷𝑐𝑜𝑠(𝜑) + 𝐸𝑠𝑖𝑛(𝜑) e 𝑗 = 𝐸𝑐𝑜𝑠(𝜑) −

𝐷𝑠𝑖𝑛(𝜑).

Pelo fato das equações polares assumirem diferentes valores em cada caso, o vetor normal

e o raio instantâneo de curvatura também serão diferentes. Substituindo a equação paramétrica

Departamento de Engenharia Mecânica

no programa “Maple” apresentado anteriormente, obtém-se o vetor normal e do raio

instantâneo para cada um dos casos apresentados

Simulação

A simulação apresentada a seguir utilizou o algoritmo apresentado em diagrama no caso da

elipse e utilizou como exemplo a equação cartesiana 2𝑋2 − 4𝑋𝑌 + 2𝑌2 − 3400𝑋 − 3800𝑌 ++5100.

clear; clc; %Coeficientes A=2; B=-4; C=2; D=-3400; E=-3800; F=5100;

if A==C psi=pi/4; else psi= atan((B/(A-

C)))*(1/2);%Angulo de Rotacao end

u = A*(cos(psi))^2 +

B*sin(psi)*cos(psi) +

C*(sin(psi))^2; v = A*(sin(psi))^2 -

B*sin(psi)*cos(psi) +

C*(cos(psi))^2; k = D*cos(psi) + E*sin(psi); j = E*cos(psi) - D*sin(psi);

V0=3; %Velocidade Inicial V(1)=V0; w0=0;%Velocidade Angular Inicial w(1)=w0; At(1)=1;%Aceleracao Tangencial

Inicial

An(1)= 0;%Aceleracao Normal

Inicial t(1)=0;%Tempo Inicial A(1)=0;%Aceleracao Inicial

tf = 10; %Tempo Final dt = .005; %Intervalo de Simulacao

%Caso u aprox. 0 if abs(u)<= 0.1 %Posicao do Vertice em

coordenadas locais X0 = -(F/k) + (j^2)/(4*v*k); Y0 = -j/(2*v); p=-k/(2*v);%2x a distancia

focal

if p>0 theta0=pi/4;%Angulo

Inicial Rho(1) =p/((1-

(cos(theta0))));%Equacao Polar else theta0=-pi/4;%Angulo

Inicial Rho(1) =

p/(1+(cos(theta0)));%Equacao Polar end theta(1)=theta0;

%Posicoes em diferentes

coordenadas Xlll(1) = Rho(1)*cos(theta0); if p>=0 Xll(1)= Xlll(1) + p/2; else Xll(1)= Xlll(1) - p/2; end Yll(1)= Rho(1)*sin(theta0); Xl(1) = Xll(1) + X0; Yl(1) = Yll(1) + Y0; X(1) = Xl(1)*cos(psi) -

Yl*sin(psi); Y(1) = Xl*sin(psi) +

Yl*cos(psi);

%Definicao do Raio de

Curvatura Instantaneo e Vetor

Normal if p>0 R(1) = RdCPUP(theta0, k,

v); N= VNPUP(theta0,k,v); else R(1) = RdCPUN(theta0, k,

v); N= VNPUN(theta0,k,v); end

%Centro Instantaneo em

Diferentes Coordenadas NXll(1) = R(1)*N(1,1)+ Xll(1); NYll(1) = R(1)*N(2,1)+ Yll(1); NXl(1) = X0 + NXll(1); NYl(1) =Y0 + NYll(1);

Departamento de Engenharia Mecânica

NX(1)= cos(psi)*NXl(1) -

sin(psi)*NYl(1); NY(1) = sin(psi)*NXl(1) +

cos(psi)*NYl(1);

%Velocidade em Diferentes

Coordenadas if p>0

Modulo0=sqrt(((0.5*k*cos(theta(1))

*sin(theta(1))/(v*((1-

cos(theta(1)))^2)))+(0.5*k*(sin(th

eta(1))/(v*(1-cos(theta(1)))))))^2

+

(((.5*k*(sin(theta(1))^2)/(v*((1-

cos(theta(1)))^2))) +

(.5*k*cos(theta(1))/(v*(1-

cos(theta(1)))))))^2); Vxll(1) =

V0*(((0.5*k*cos(theta(1))*sin(thet

a(1))/(v*((1-

cos(theta(1)))^2)))+(0.5*k*(sin(th

eta(1))/(v*(1-

cos(theta(1))))))))/Modulo0; Vyll(1) =

V0*(((.5*k*(sin(theta(1))^2)/(v*((

1-cos(theta(1)))^2))) +

(.5*k*cos(theta(1))/(v*(1-

cos(theta(1)))))))/Modulo0; else Modulo0=sqrt(((-

0.5*k*cos(theta(1))*sin(theta(1))/

(v*((1+cos(theta(1)))^2)))+(0.5*k(

sin(theta(1))/(v*(1+cos(theta(1)))

))))^2 + (((-

.5*k*(sin(theta(1))^2)/(v*((1+cos(

theta(1)))^2))) -

(.5*k*cos(theta(1))/(v*(1+cos(thet

a(1)))))))^2); Vxll(1) = V0*(((-

0.5*k*cos(theta(1))*sin(theta(1))/

(v*((1+cos(theta(1)))^2)))+(0.5*k(

sin(theta(1))/(v*(1+cos(theta(1)))

)))))/Modulo0; Vyll(1) = V0*(((-

.5*k*(sin(theta(1))^2)/(v*((1+cos(

theta(1)))^2))) -

(.5*k*cos(theta(1))/(v*(1+cos(thet

a(1)))))))/Modulo0; end Vx(1)=cos(psi)*Vxll(1) -

sin(psi)*Vyll(1); Vy(1)=sin(psi)*Vxll(1) +

cos(psi)*Vyll(1);

for i=1:1:tf/dt t(i+1)=i*dt; if t(i+1)<10 At(i+1) = 1; elseif t(i+1)>=10 && t(i+1)<30

At(i+1) = 2; else At(i+1) =0; end

V(i+1)=V(i)+At(i)*dt; An(i+1)=(V(i))^2/R(i); A(i+1)=sqrt((An(i+1))^2 +

(At(i+1)^2)); w(i+1)=V(i)/R(i); theta(i+1)=theta(i)+w(i)*dt +

At(i+1)/R(i)*dt^2; if p>0 R(i+1) =

RdCPUP(theta(i+1), k, v); Rho(i+1) = p/(1-

(cos(theta(i+1)))); else R(i+1) =

RdCPUN(theta(i+1), k, v); Rho(i+1) = -

p/(1+(cos(theta(i+1)))); end

Xlll(i+1) =

Rho(i+1)*cos(theta(i+1)); if p>0 Xll(i+1)= Xlll(i+1) + p/2; else Xll(i+1)= Xlll(i+1) - p/2; end Yll(i+1) =

Rho(i+1)*sin(theta(i+1)); Xl(i+1) = X0 + Xll(i+1); Yl(i+1) = Y0 + Yll(i+1); X(i+1)= cos(psi)*Xl(i+1) -

sin(psi)*Yl(i+1); Y(i+1) = sin(psi)*Xl(i+1) +

cos(psi)*Yl(i+1);

if p>0 N= VNPUP(theta(i+1),k,v); else N= VNPUN(theta(i+1),k,v); end

NXll(i+1) = R(i+1)*N(1,1)+

Xll(i+1); NYll(i+1) = R(i+1)*N(2,1)+

Yll(i+1); NXl(i+1) = X0 + NXll(i+1); NYl(i+1) =Y0 + NYll(i+1); NX(i+1)= cos(psi)*NXl(i+1) -

sin(psi)*NYl(i+1); NY(i+1) = sin(psi)*NXl(i+1) +

cos(psi)*NYl(i+1);

if p>0

ModuloI=sqrt(((0.5*k*cos(theta(i+1

Departamento de Engenharia Mecânica

))*sin(theta(i+1))/(v*((1-

cos(theta(i+1)))^2)))+(0.5*k*(sin(

theta(i+1))/(v*(1-

cos(theta(i+1)))))))^2 +

(((.5*k*(sin(theta(i+1))^2)/(v*((1

-cos(theta(i+1)))^2))) +

(.5*k*cos(theta(i+1))/(v*(1-

cos(theta(i+1)))))))^2); Vxll(i+1) =

V(i+1)*(((0.5*k*cos(theta(i+1))*si

n(theta(i+1))/(v*((1-

cos(theta(i+1)))^2)))+(0.5*k*(sin(

theta(i+1))/(v*(1-

cos(theta(i+1))))))))/ModuloI; Vyll(i+1) =

V(i+1)*(((.5*k*(sin(theta(i+1))^2)

/(v*((1-cos(theta(i+1)))^2))) +

(.5*k*cos(theta(i+1))/(v*(1-

cos(theta(i+1)))))))/ModuloI; else ModuloI=sqrt(((-

0.5*k*cos(theta(i+1))*sin(theta(i+

1))/(v*((1+cos(theta(i+1)))^2)))+(

0.5*k(sin(theta(i+1))/(v*(1+cos(th

eta(i+1)))))))^2 + (((-

.5*k*(sin(theta(i+1))^2)/(v*((1+co

s(theta(i+1)))^2))) -

(.5*k*cos(theta(i+1))/(v*(1+cos(th

eta(i+1)))))))^2); Vxll(i+1) = V(i+1)*(((-

0.5*k*cos(theta(i+1))*sin(theta(i+

1))/(v*((1+cos(theta(i+1)))^2)))+(

0.5*k(sin(theta(i+1))/(v*(1+cos(th

eta(i+1))))))))/ModuloI; Vyll(i+1) = V(i+1)*(((-

.5*k*(sin(theta(i+1))^2)/(v*((1+co

s(theta(i+1)))^2))) -

(.5*k*cos(theta(i+1))/(v*(1+cos(th

eta(i+1)))))))/ModuloI; end Vx(i+1)=cos(psi)*Vxll(i+1) -

sin(psi)*Vyll(i+1); Vy(i+1)=sin(psi)*Vxll(i+1) +

cos(psi)*Vyll(i+1);

end

%Caso v aprox. 0 else %Posicao do Vertice em

coordenadas locais X0 = -k/(2*u); Y0 = -(F/j) + (k^2)/(4*u*j); p=-j/(2*u);%2x a distancia

focal

if p>0 theta0=-.7*pi;%Angulo

Inicial Rho(1) = p/(1-

(sin(theta0)));%Equacao Polar

else theta0=-pi/4;%Angulo

Inicial Rho(1) =

p/(1+(sin(theta0)));%Equacao Polar end theta(1)=theta0;

%Posicao em diferentes

coordenadas Ylll(1) = Rho(1)*sin(theta0); if p>=0 Yll(1)= Ylll(1) + p/2; else Yll(1)=Ylll(1) - p/2; end

Xll(1) = Rho(1)*cos(theta0); Xl(1) = Xll(1) + X0; Yl(1) = Yll(1) + Y0; X(1) = Xl(1)*cos(psi) -

Yl*sin(psi); Y(1) = Xl*sin(psi) +

Yl*cos(psi);

%Definicao do Raio de

Curvatura Instantaneo e Vetor

Normal if p>0 R(1) = RdCPVP(theta0, j,

u); N= VNPVP(theta0,j,u); else R(1) = RdCPVN(theta0, j,

u); N= VNPVN(theta0,j,u); end %Centro Instantaneo em

Diferentes Coordenadas NXll(1) = R(1)*N(1,1)+ Xll(1); NYll(1) = R(1)*N(2,1)+ Yll(1); NXl(1) = X0 + NXll(1); NYl(1) =Y0 + NYll(1); NX(1)= cos(psi)*NXl(1) -

sin(psi)*NYl(1); NY(1) = sin(psi)*NXl(1) +

cos(psi)*NYl(1);

%Velocidade em diferentes

coordenadas if p>0 Modulo0=sqrt(((-

0.5*j*cos(theta(1))*sin(theta(1))/

(u*((1-sin(theta(1)))^2)))-

(0.5*j*(cos(theta(1))/(u*(1-

sin(theta(1)))))))^2 +

(((.5*j*(cos(theta(1))^2)/(u*((1-

sin(theta(1)))^2))) +

(.5*j*sin(theta(1))/(u*(1-

sin(theta(1))))))^2));

Departamento de Engenharia Mecânica

Vxll(1) =

V0*(((.5*j*(cos(theta(1))^2)/(u*((

1-sin(theta(1)))^2))) +

(.5*j*sin(theta(1))/(u*(1-

sin(theta(1)))))))/Modulo0; Vyll(1) = V0*(((-

0.5*j*cos(theta(1))*sin(theta(1))/

(u*((1-sin(theta(1)))^2)))-

(0.5*j*(cos(theta(1))/(u*(1-

sin(theta(1))))))))/Modulo0; else

Modulo0=sqrt(((0.5*j*cos(theta(1))

*sin(theta(1))/(u*((1+sin(theta(1)

))^2)))-

(0.5*j*(cos(theta(1))/(u*(1+sin(th

eta(1)))))))^2 +

(((.5*j*(cos(theta(1))^2)/(u*((1+s

in(theta(1)))^2))) +

(.5*j*sin(theta(1))/(u*(1+sin(thet

a(1))))))^2)); Vxll(1) =

V0*(((.5*j*(cos(theta(1))^2)/(u*((

1+sin(theta(1)))^2))) +

(.5*j*sin(theta(1))/(u*(1+sin(thet

a(1)))))))/Modulo0; Vyll(1) =

V0*(((0.5*j*cos(theta(1))*sin(thet

a(1))/(u*((1+sin(theta(1)))^2)))-

(0.5*j*(cos(theta(1))/(u*(1+sin(th

eta(1))))))))/Modulo0; end

Vx(1)=cos(psi)*Vxll(1) -

sin(psi)*Vyll(1); Vy(1)=sin(psi)*Vxll(1) +

cos(psi)*Vyll(1);

for i=1:1:tf/dt t(i+1)=i*dt; if t(i+1)<10 At(i+1) = 5; elseif t(i+1)>=10 && t(i+1)<30 At(i+1) = 4; else At(i+1) = 2; end

V(i+1)=V(i)+At(i)*dt; An(i+1)=(V(i))^2/R(i); A(i+1)=sqrt((An(i+1))^2 +

(At(i+1)^2)); w(i+1)=V(i)/R(i); theta(i+1)=theta(i)+w(i)*dt +

At(i+1)/R(i)*dt^2;

if p>0 R(i+1) = RdCPVP(theta(i+1),

j, u);

Rho(i+1) = p/(1-

(sin(theta(i+1)))); else R(i+1) =

RdCPVN(theta(i+1), j, u); Rho(i+1) = -

p/(1+(sin(theta(i+1)))); end

Ylll(i+1) =

Rho(i+1)*sin(theta(i+1)); if p>=0 Yll(i+1)= Ylll(i+1) +p/2 ; else Yll(i+1)=Ylll(i+1) - p/2; end

Xll(i+1) =

Rho(i+1)*cos(theta(i+1)); Xl(i+1) = X0 + Xll(i+1); Yl(i+1) = Y0 + Yll(i+1); X(i+1)= cos(psi)*Xl(i+1) -

sin(psi)*Yl(i+1); Y(i+1) = sin(psi)*Xl(i+1) +

cos(psi)*Yl(i+1);

if p>0 N= VNPVP(theta(i+1),j,u); else N= VNPVN(theta(i+1),j,u); end

NXll(i+1) = R(i+1)*N(1,1)+

Xll(i+1); NYll(i+1) = R(i+1)*N(2,1)+

Yll(i+1); NXl(i+1) = X0 + NXll(i+1); NYl(i+1) =Y0 + NYll(i+1); NX(i+1)= cos(psi)*NXl(i+1) -

sin(psi)*NYl(i+1); NY(i+1) = sin(psi)*NXl(i+1) +

cos(psi)*NYl(i+1);

if p>0 ModuloI=sqrt(((-

0.5*j*cos(theta(i+1))*sin(theta(i+

1))/(u*((1-sin(theta(i+1)))^2)))-

(0.5*j*(cos(theta(i+1))/(u*(1-

sin(theta(i+1)))))))^2 +

(((.5*j*(cos(theta(i+1))^2)/(u*((1

-sin(theta(i+1)))^2))) +

(.5*j*sin(theta(i+1))/(u*(1-

sin(theta(i+1))))))^2)); Vxll(i+1) =

V(i+1)*(((.5*j*(cos(theta(i+1))^2)

/(u*((1-sin(theta(i+1)))^2))) +

(.5*j*sin(theta(i+1))/(u*(1-

sin(theta(i+1)))))))/ModuloI; Vyll(i+1) = V(i+1)*(((-

0.5*j*cos(theta(i+1))*sin(theta(i+

Departamento de Engenharia Mecânica

1))/(u*((1-sin(theta(i+1)))^2)))-

(0.5*j*(cos(theta(i+1))/(u*(1-

sin(theta(i+1))))))))/ModuloI; else

ModuloI=sqrt(((0.5*j*cos(theta(i+1

))*sin(theta(i+1))/(u*((1+sin(thet

a(i+1)))^2)))-

(0.5*j*(cos(theta(i+1))/(u*(1+sin(

theta(i+1)))))))^2 +

(((.5*j*(cos(theta(i+1))^2)/(u*((1

+sin(theta(i+1)))^2))) +

(.5*j*sin(theta(i+1))/(u*(1+sin(th

eta(i+1))))))^2)); Vxll(i+1) =

V(i+1)*(((.5*j*(cos(theta(i+1))^2)

/(u*((1+sin(theta(i+1)))^2))) +

(.5*j*sin(theta(i+1))/(u*(1+sin(th

eta(i+1)))))))/ModuloI; Vyll(i+1) =

V(i+1)*(((0.5*j*cos(theta(i+1))*si

n(theta(i+1))/(u*((1+sin(theta(i+1

)))^2)))-

(0.5*j*(cos(theta(i+1))/(u*(1+sin(

theta(i+1))))))))/ModuloI; end

Vx(i+1)=cos(psi)*Vxll(i+1) -

sin(psi)*Vyll(i+1); Vy(i+1)=sin(psi)*Vxll(i+1) +

cos(psi)*Vyll(i+1); end

end

Departamento de Engenharia Mecânica

4) Hipérbole

A Hipérbole é uma curva que possui duas partes espelhadas pelo seu eixo de simetria.

Devido esse fato, para que se pudesse parametrizar essa cônica, foi necessário a escolha de uma

de suas metades para poder descrever a trajetória da partícula orientada. O modelo da hipérbole

é bem semelhante ao da elipse, utilizando a mesma ordem de transformações de coordenadas

e a mesma equação para determinação do centro. Com isso, o modelo da hipérbole é descrito

por dois casos:

Caso 1: 𝑋2

𝑎2 −𝑌2

𝑏2 = 1.

O eixo de simetria para esse caso é paralelo ao eixo Y’’ e a metade escolhida para a

parametrização foi à esquerda, embora não haja problema quanto a escolha da outra parte. A

equação polar obtida é dada por 𝜌(𝜃) =𝑏2

𝑎+𝑐∗𝑐𝑜𝑠(𝜃) , onde 𝑐 = √𝑎2 + 𝑏2. 𝑎 e 𝑏 são os mesmos

da elipse, salvo o sinal de menos, sendo substituído pelo operador módulo, já que a combinação

dessas variáveis com 𝑢 e 𝑣 faz com que haja os dois casos analisados.

Caso 2: −𝑋2

𝑎2 +𝑌2

𝑏2 = 1.

O eixo de simetria para esse caso é paralelo ao eixo X’’ e a metade escolhida para a

parametrização foi a inferior, embora não haja problema quanto a escolha da outra parte. A

equação polar obtida é dada por 𝜌(𝜃) =𝑏2

𝑎+𝑐∗𝑠𝑖𝑛(𝜃). As outras variáveis utilizadas são

semelhantes ao caso 1 e ao caso da elipse.

Simulação

A simulação apresentada a seguir utilizou o algoritmo apresentado em diagrama no caso da

elipse e utilizou como exemplo a equação cartesiana 2𝑋2 + 9𝑋𝑌 + 2𝑌2 − 4000𝑋 − 3000𝑌 ++7000.

clear; clc; %Coeficientes A=2; B=9; C=2; D=-4000; E=-3000; F=7000;

%Variaveis Iniciais V0=3; %Velocidade Inicial em m/s X0= -((-B * E + 2 * C * D + B * D)

/ A / (2 * C + B)) /

0.2e1;%Posição do centro no eixo X Y0= -E / (2 * C + B);%Posição do

centro no eixo Y Fl = A*X0^2 + B*X0*Y0 + C*Y0^2 +

D*X0 + E*Y0 + F;

if A==C psi=pi/4; else psi= atan((B/(A-

C)))*(1/2);%Angulo de Rotacao end

u = A*(cos(psi))^2 +

B*sin(psi)*cos(psi) +

C*(sin(psi))^2; v = A*(sin(psi))^2 -

B*sin(psi)*cos(psi) +

C*(cos(psi))^2; a= sqrt(abs(Fl/u));%Coeficiente de

Xll b= sqrt(abs(Fl/v));%Coeficiente de

Yll

Departamento de Engenharia Mecânica

c=sqrt(a^2 + b^2);

%Caso (X^2/a^2)-(Y^2/b^2)=1 if (u>0 && v<0 && Fl<0) || (u<0 &&

v>0 && Fl>0) theta0= -pi/4;%Angulo Inicial theta(1)=theta0; w0=0;%Velocidade Angular

Inicial w(1)=w0; At(1)=1;%Aceleracao Tangencial

Inicial

Rho(1)=(b^2)/(a+c*cos(theta(1)));

%Posicoes em direfentes

coordenadas Xlll(1)=Rho(1)*cos(theta(1)); Xll(1)= Xlll(1)-c; Yll(1)= Rho(1)*sin(theta(1)); Xl(1) = cos(psi)*Xll(1) -

sin(psi)*Yll(1); Yl(1) = sin(psi)*Xll(1) +

cos(psi)*Yll(1); X(1)= X0 + Xl(1); Y(1) = Y0 + Yl(1);

R(1) = RdCHP(theta(1), a,

b,c);%Raio de Curvatura

Instantaneo N= VNHP(theta(1),a,b,c);%Vetor

Normal NXll(1) = R(1)*N(1,1)+ Xll(1); NYll(1) = R(1)*N(2,1)+ Yll(1); NXl(1) = cos(psi)*NXll(1) -

sin(psi)*NYll(1); NYl(1) = sin(psi)*NXll(1) +

cos(psi)*NYll(1); NX(1)= X0 + NXl(1); NY(1) = Y0 + NYl(1);

%Velocidades em Diferentes

Coordenadas Modulo0=(a + c *

cos(theta(1))) *

sin(theta(1))/sqrt((b ^ 2 / (a + c

* cos(theta(1))) ^ 2 *

cos(theta(1)) * c * sin(theta(1))

- b ^ 2 / (a + c * cos(theta(1)))

* sin(theta(1)))^2 + (b ^ 2 / (a +

c * cos(theta(1))) ^ 2 *

sin(theta(1)) ^ 2 * c + b ^ 2 / (a

+ c * cos(theta(1))) *

cos(theta(1)))^2); Vxll(1) = V0*(b ^ 2 / (a + c *

cos(theta(1))) ^ 2 * cos(theta(1))

* c * sin(theta(1)) - b ^

2)/Modulo0;

Vyll(1) = V0*(b ^ 2 / (a + c *

cos(theta(1))) ^ 2 * sin(theta(1))

^ 2 * c + b ^ 2)/Modulo0; Vx(1)=cos(psi)*Vxll(1) -

sin(psi)*Vyll(1); Vy(1)=sin(psi)*Vxll(1) +

cos(psi)*Vyll(1); V(1)=sqrt(Vx^2 + Vy^2);

An(1)= 0;%Aceleracao Normal

Inicial t(1)=0;%Tempo Inicial A(1)=0;%Modulo da Aceleracao

Inicial

tf = 60; %Tempo Final dt = .005; %Intervalo de

Simulacao

for i=1:1:tf/dt t(i+1)=i*dt; if t(i+1)<5 At(i+1) = 1; elseif t(i+1)>=5 &&

t(i+1)<7 At(i+1) = .5; else At(i+1) = 0; end

V(i+1)= V(i)+At(i)*dt; An(i+1)=(V(i))^2/R(i); A(i+1)= sqrt((An(i+1))^2 +

(At(i+1)^2)); w(i+1)= V(i)/R(i); theta(i+1)=

theta(i)+w(i)*dt +

At(i+1)/R(i)*dt^2; Rho(i+1)=

(b^2)/(a+c*cos(theta(i+1))); R(i+1) = RdCHP(theta(i+1),

a, b,c);

Xlll(i+1)=Rho(i+1)*cos(theta(i+1))

; Xll(i+1)= Xlll(i+1)-c; Yll(i+1)=

Rho(i+1)*sin(theta(i+1)); Xl(i+1) =

cos(psi)*Xll(i+1) -

sin(psi)*Yll(i+1); Yl(i+1) =

sin(psi)*Xll(i+1) +

cos(psi)*Yll(i+1); X(i+1)= X0 + Xl(i+1); Y(i+1) = Y0 + Yl(i+1);

N = VNHP(theta(i+1), a,

b,c);

Departamento de Engenharia Mecânica

NXll(i+1) = R(i+1)*N(1,1)+

Xll(i+1); NYll(i+1) = R(i+1)*N(2,1)+

Yll(i+1); NXl(i+1) =

cos(psi)*NXll(i+1) -

sin(psi)*NYll(i+1); NYl(i+1) =

sin(psi)*NXll(i+1) +

cos(psi)*NYll(i+1); NX(i+1)= X0 + NXl(i+1); NY(i+1) = Y0 + NYl(i+1);

ModuloI= (a + c *

cos(theta(i+1))) *

sin(theta(i+1))/ sqrt((b ^ 2 / (a

+ c * cos(theta(i+1))) ^ 2 *

cos(theta(i+1)) * c *

sin(theta(i+1)) - b ^ 2 / (a + c *

cos(theta(i+1))) *

sin(theta(i+1)))^2 + (b ^ 2 / (a +

c * cos(theta(i+1))) ^ 2 *

sin(theta(i+1)) ^ 2 * c + b ^ 2 /

(a + c * cos(theta(i+1))) *

cos(theta(i+1)))^2); Vxll(i+1) = V(i+1)*(b ^ 2

/ (a + c * cos(theta(i+1))) ^ 2 *

cos(theta(i+1)) * c *

sin(theta(i+1)) - b ^ 2)/ModuloI; Vyll(i+1) = V(i+1)*(b ^ 2

/ (a + c * cos(theta(i+1))) ^ 2 *

sin(theta(i+1)) ^ 2 * c + b ^

2)/ModuloI; Vx(i+1)=cos(psi)*Vxll(1) -

sin(psi)*Vyll(1); Vy(i+1)=sin(psi)*Vxll(1) +

cos(psi)*Vyll(1);

end

%Caso -(X^2/a^2)+(Y^2/b^2)=1 else theta0= 0; %Angulo Inicial theta(1)=theta0; w0=0;%Velocidade Angular

Inicial w(1)=w0; At(1)=1;%Aceleracao Tangencial

Inicial

Rho(1)=b^2/(a+c*sin(theta(1)));

%Posicoes em direfentes

coordenadas Xll(1)=Rho(1)*cos(theta(1)); Ylll(1)= Rho(1)*sin(theta(1)); Yll(1)= Ylll(1)-c; Xl(1) = cos(psi)*Xll(1) -

sin(psi)*Yll(1);

Yl(1) = sin(psi)*Xll(1) +

cos(psi)*Yll(1); X(1)= X0 + Xl(1); Y(1) = Y0 + Yl(1);

R(1) = RdCHS(theta(1), a,

b,c);%Raio de Curvatura

Instantaneo N= VNHS(theta(1),a,b,c);%Vetor

Normal NXll(1) = R(1)*N(1,1)+ Xll(1); NYll(1) = R(1)*N(2,1) +

Yll(1); NXl(1) = cos(psi)*NXll(1) -

sin(psi)*NYll(1); NYl(1) = sin(psi)*NXll(1) +

cos(psi)*NYll(1); NX(1)= X0 + NXl(1); NY(1) = Y0 + NYl(1);

%Velocidades em Diferentes

Coordenadas Modulo0=sqrt((-

(b^2)*(cos(theta(1))^2)*c/((a+c*si

n(theta(1)))^2)-

(b^2*sin(theta(1))/(a+c*sin(theta(

1)))))^2 + (((-

(b^2)*sin(theta(1))*c*cos(theta(1)

))/((a+c*sin(theta(1)))^2)) +

b^2*cos(theta(1))/(a+c*sin(theta(1

))))^2); Vxll(1) = V0*(-

(b^2)*(cos(theta(1))^2)*c/((a+c*si

n(theta(1)))^2)-

(b^2*sin(theta(1))/(a+c*sin(theta(

1)))))/Modulo0; Vyll(1) = V0*(-

(b^2)*sin(theta(1))*c*cos(theta(1)

)/((a+c*sin(theta(1)))^2) +

b^2*cos(theta(1))/(a+c*sin(theta(1

))))/Modulo0; Vx(1)=cos(psi)*Vxll(1) -

sin(psi)*Vyll(1); Vy(1)=sin(psi)*Vxll(1) +

cos(psi)*Vyll(1);

An(1)= 0;%Aceleracao Normal

Inicial t(1)=0;%Tempo Inicial A(1)=0;%Modulo da Aceleracao

Inicial

tf = 60; %Tempo Final dt = .005; %Intervalo de

Simulacao

for i=1:1:tf/dt t(i+1)=i*dt; if t(i+1)<5 At(i+1) = 1;

Departamento de Engenharia Mecânica

elseif t(i+1)>=5 &&

t(i+1)<7 At(i+1) = .5; else At(i+1) = 0; end

V(i+1)= V(i)+At(i)*dt; An(i+1)=(V(i))^2/R(i); A(i+1)= sqrt((An(i+1))^2 +

(At(i+1)^2)); w(i+1)= V(i)/R(i); theta(i+1)=

theta(i)+w(i)*dt +

At(i+1)/R(i)*dt^2; Rho(i+1)=

b^2/(a+c*sin(theta(i+1))); R(i+1) = RdCHS(theta(i+1),

a, b,c);

Xll(i+1)=Rho(i+1)*cos(theta(i+1)); Ylll(i+1)=

Rho(i+1)*sin(theta(i+1)); Yll(i+1)= Ylll(i+1) -c; Xl(i+1) =

cos(psi)*Xll(i+1) -

sin(psi)*Yll(i+1); Yl(i+1) =

sin(psi)*Xll(i+1) +

cos(psi)*Yll(i+1); X(i+1)= X0 + Xl(i+1); Y(i+1) = Y0 + Yl(i+1);

N = VNHS(theta(i+1), a,

b,c); NXll(i+1) = R(i+1)*N(1,1)+

Xll(i+1);

NYll(i+1) =

R(i+1)*N(2,1)+Yll(i+1); NXl(i+1) =

cos(psi)*NXll(i+1) -

sin(psi)*NYll(i+1); NYl(i+1) =

sin(psi)*NXll(i+1) +

cos(psi)*NYll(i+1); NX(i+1)= X0 + NXl(i+1); NY(i+1) = Y0 + NYl(i+1);

ModuloI=sqrt((-

(b^2)*(cos(theta(i+1))^2)*c/((a+c*

sin(theta(i+1)))^2)-

(b^2*sin(theta(i+1))/(a+c*sin(thet

a(i+1)))))^2 + (((-

(b^2)*sin(theta(i+1))*c*cos(theta(

i+1)))/((a+c*sin(theta(i+1)))^2))

+

b^2*cos(theta(i+1))/(a+c*sin(theta

(i+1))))^2); Vxll(i+1) =V(1+1)*(-

(b^2)*(cos(theta(i+1))^2)*c/((a+c*

sin(theta(i+1)))^2)-

(b^2*sin(theta(i+1))/(a+c*sin(thet

a(i+1)))))/ModuloI; Vyll(i+1) = V(i+1)*(-

(b^2)*sin(theta(i+1))*c*cos(theta(

i+1))/((a+c*sin(theta(i+1)))^2) +

b^2*cos(theta(i+1))/(a+c*sin(theta

(i+1))))/ModuloI; Vx(i+1)=cos(psi)*Vxll(1) -

sin(psi)*Vyll(1); Vy(i+1)=sin(psi)*Vxll(1) +

cos(psi)*Vyll(1); end end

Departamento de Engenharia Mecânica

Conclusão

A obtenção de todas as informações necessárias para o futuro estudo das trajetórias

ótimas aproximadas por trechos de curvas cônicas foram encontradas com sucesso. Além disso,

o programa foi testado para diferentes equações paramétricas, obtendo resultados esperados.

Diante do sucesso dessa primeira etapa, tem-se como objetivo o avanço do projeto na

elaboração de trajetórias ótimas a partir do material desenvolvido nesse trabalho.

Bibliografia

1. VENTURI, Jacir J., Cônicas e Quádricas 5.a edição, páginas 23-135.

2. GINSBERG, Jerry, Engineering Dynamics, Cambridge, cap 2.1 – 2.1.2.

3. da CRUZ, Luiz Francisco, Cálculo Vetorial e Geometria Analítica, Departamento de

Matemática – Unesp/Bauru, cap 9.