DM de potenciais descontínuos: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
 
(15 revisões intermediárias por 2 usuários não estão sendo mostradas)
Linha 4: Linha 4:
[[File:fluxogram.png|thumb|Fluxograma de um programa simples usando o método de evento dirigido com a otimização aqui explicitada.]]
[[File:fluxogram.png|thumb|Fluxograma de um programa simples usando o método de evento dirigido com a otimização aqui explicitada.]]
===Determinação do tempo de colisão===
===Determinação do tempo de colisão===
Os objetos a serem usados para o cálculo do tempo de colisão entre um par de partículas <math>i, j</math> serão discos de raio <math>\sigma_i</math>, <math>\sigma_j</math>, de distância denotada por <math>\sigma</math>. Portanto, segue que a condição de colisão é:  
Os objetos a serem usados para o cálculo do tempo de colisão entre um par de partículas <math>i, j</math> serão discos de raio <math>\sigma_i</math>, <math>\sigma_j</math>, cuja soma denotamos por <math>\sigma</math>. Portanto, segue que a condição de colisão é:  
:<math>|\vec{r_i}(t + dt_{ij}) - \vec{r_j}(t + dt_{ij})| = \sigma</math>
:<math>|\vec{r_i}(t + dt_{ij}) - \vec{r_j}(t + dt_{ij})| = \sigma</math>
Com <math>r_i</math> sendo o vetor posição da partícula <math>i</math> e <math>dt_{ij}</math> o tempo de colisão entre as partículas <math>i, j</math>. Tal condição nos leva a determinação de <math>dt_{ij}</math> a partir da expressão:
Com <math>\vec{r_i}</math> sendo o vetor posição da partícula <math>i</math> e <math>dt_{ij}</math> o tempo de colisão entre as partículas <math>i, j</math>. Tal condição nos leva a determinação de <math>dt_{ij}</math> a partir da expressão:
:<math>
:<math>
dt_{ij} =
dt_{ij} =
Linha 18: Linha 18:


===Mudança de velocidade em uma colisão elástica===
===Mudança de velocidade em uma colisão elástica===
Para fazer a mudança de velocidades temos que considerar o caso de colisão elástica entre as partículas <math> i, j </math>, sendo impulso dado por:
Para fazer a mudança de velocidades temos que considerar o caso de colisão elástica entre as partículas <math> I, J </math>, sendo impulso dado por:
:<math> \vec{J} = \frac{2m_im_j(\Delta \vec{r} . \Delta \vec{v})}{\sigma^2(m_i + m_j)}\Delta \vec{r} </math>.
:<math> \vec{J} = \frac{2m_im_j(\Delta \vec{r} . \Delta \vec{v})}{\sigma^2(m_i + m_j)}\Delta \vec{r} </math>.
Assim, a variação de velocidades pode ser determinada por:
Assim, a variação de velocidades pode ser determinada por:
Linha 27: Linha 27:


===Implementação computacional===
===Implementação computacional===
Segue a implementação computacional, na linguagem de programação C, a função utilizada para o cálculo de <math>dt_{ij}</math>.
Segue a implementação computacional, na linguagem de programação C, a função utilizada para o cálculo de <math>dt_{ij}</math>. Importante notar que na formulação da funções-exemplo foram utilizadas condições de contorno periódicas no sistema.
<pre>
<pre>
void calc_dt(double *xx, double *yy, double *vx, double *vy, double *deltat, int *colide, int i, int j){
void calc_dt_ij(double *xx, double *yy, double *vx, double *vy, double *dt_ij, int *colide, int i, int j){


   double delta_t = INF;
//dt_ij é um vetor que armazena os menores valor de dt_ij para cada i
//colide é um vetor que armazena o valor de j correspondente a partícula que i colide no menor tempo dt_ij
 
   double delta_t_ij = INF;
   double dx, dy, dvx, dvy;
   double dx, dy, dvx, dvy;
   double drdr, dvdv, drdv;
   double drdr, dvdv, drdv;
Linha 53: Linha 56:


   if(d > 0 && drdv < 0){
   if(d > 0 && drdv < 0){
     delta_t = (-1)*(drdv + sqrt(d))/(dvdv);
     delta_t_ij = (-1)*(drdv + sqrt(d))/(dvdv);
       if(delta_t < deltat[i]){
       if(delta_t_ij < dt_ij[i]){
deltat[i] = delta_t;
dt_ij[i] = delta_t_ij;
colide[i] = j;
  colide[i] = j;
       }
       }
   }
   }
Linha 91: Linha 94:
   deltavy = Jy[i];
   deltavy = Jy[i];


   //TROCA VELOCIDADES DA PARTICULA i
   //TROCA VELOCIDADES DA PARTICULA I
   vx[i] -= deltavx;
   vx[i] -= deltavx;
   vy[i] -= deltavy;
   vy[i] -= deltavy;


   //TROCA VELOCIDADES DA PARTICULA j
   //TROCA VELOCIDADES DA PARTICULA J
   vx[j] += deltavx;
   vx[j] += deltavx;
   vy[j] += deltavy;
   vy[j] += deltavy;
Linha 103: Linha 106:
Após alterar as velocidades basta avançar o sistema em <math>dt_{min}</math>. Segue um exemplo de função.
Após alterar as velocidades basta avançar o sistema em <math>dt_{min}</math>. Segue um exemplo de função.
<pre>
<pre>
void pos(double *xx, double *yy, double *vx, double *vy, double delta_t){
void pos(double *xx, double *yy, double *vx, double *vy, double dt_min){


   int i;
   int i;


   for(i = 0; i < (NP); i++){
   for(i = 0; i < (NP); i++){
     xx[i] = xx[i] + vx[i]*delta_t;
     xx[i] = xx[i] + vx[i]*dt_min;
     yy[i] = yy[i] + vy[i]*delta_t;
     yy[i] = yy[i] + vy[i]*dt_min;
   }
   }


Linha 115: Linha 118:
</pre>
</pre>
[[File:pot_descontinuo.gif|400px|thumb|Plot exemplo de resultado de simulação. O par de partículas destacadas em vermelho são as colisoras.]]
[[File:pot_descontinuo.gif|400px|thumb|Plot exemplo de resultado de simulação. O par de partículas destacadas em vermelho são as colisoras.]]
Seguindo o fluxograma apresentado e utilizando as funções disponíveis chega-se, por exemplo, na animação ao lado, que mostra o decorrer da simulação conforme ocorrem as colisões. <br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>
Seguindo o fluxograma apresentado e utilizando as funções disponíveis chega-se, por exemplo, na animação ao lado, que mostra o decorrer da simulação conforme ocorrem as colisões. <br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>


==Adição do campo gravitacional==
==Adição do campo gravitacional==
A adição de um campo gravitacional uniforme nesta simulação é bastante simples, basta definir uma nova variável, a qual aqui será chamada de g, e atribuir um valor de sua escolha a ela. Vale lembrar que a força gravitacional aqui computada gera um campo uniforme, ou seja, é radial em relação a cada partícula do sistema, de modo que só há força atuando na direção <math>y</math>, no sentido negativo de <math>y</math>.
A adição de um campo gravitacional uniforme nesta simulação é bastante simples, basta definir uma nova variável, a qual aqui será chamada de g, e atribuir um valor de sua escolha a ela.<br>
Vale lembrar que campo gravitacional aqui computado gera uma força central, ou seja, é radial em relação a cada partícula do sistema, de modo que só há força atuando na direção <math>y</math>, no sentido negativo de <math>y</math>.<br>
Também é importante lembrar que aqui é desconsiderado a interação gravitacional entre quaisquer duas partículas, uma vez que assume-se que massas são pequenas e desprezíveis em frente a massa que forma o campo em que estão imersas.<br>
===Cálculo do tempo de colisão===
Uma preocupação que alguém pode ter com a adição de um campo gravitacional na simulação seria de que seria necessário mudar o cálculo do tempo mínimo entre as colisões. No entanto, como será mostrado a seguir, isso não é necessário e o cálculo permanece o mesmo.<br>
De maneira geral, sem cálculos, o que deve ser feito é abrir cada vetor de deslocamento, adicionando um termo de tempo quadrádico ligado à aceleração gravitacional do campo escolhido, e elevar ao quadrado a condição de que ocorra a colisão. A partir daí tudo que se é necessário fazer é trabalhar a álgebra, que no fim, eliminará todos os termos com a aceleração. Para verificar o cálculo completo veja na página [[Cálculo do tempo de colisão com aceleração]].<br>


Como este campo só atua em <math>y</math>, basta modificar as linha de código acima para adicionar um termo de tempo quadrático à variação da posição assim a linha de código fica:
=== Implementação ===
Como este campo só atua na componente <math>y</math>, basta modificar as linha de código da seção anterior para adicionar um termo de tempo quadrático à variação da posição assim a linha de código fica:
<pre>
<pre>
yy[i] = yy[i] + vy[i]*delta_t - g*delta_t*delta_t/2;
yy[i] = yy[i] + vy[i]*delta_t - g*delta_t*delta_t/2;
Linha 129: Linha 138:
</pre>
</pre>
Mais abaixo é possível ver uma animação de como fica o sistema com essas mudanças feitas.
Mais abaixo é possível ver uma animação de como fica o sistema com essas mudanças feitas.
[[File:pot_desc_gravidade.gif|400px|thumb|Plot exemplo de resultado de simulação com um campo gravitacional uniforme de -2. O par de partículas destacadas em vermelho são as que colidiram.]]
[[File:pot_desc_gravidade.gif|400px|thumb|Plot exemplo de resultado de simulação com um campo gravitacional uniforme de g=2. O par de partículas destacadas em vermelho são as que colidiram.]]

Edição atual tal como às 21h44min de 12 de outubro de 2016

Dinâmica molecular de potenciais descontínuos é uma abordagem computacional utilizada para determinar o movimento de partículas duras que só interagem por forças de contato. Assim, fica evidente a diferença entre o potencial Lennard-Jones pois este se baseia em uma interação de curto alcance, como é mostrado em DM: um primeiro programa. Para entender como as colisões ocorrem, conhecer a forma do potencial a ser estudado é vital. Como estamos considerando corpos rígidos, ou seja, que não sofrem deformação, percebe-se que a força de contato entre as partículas será infinita e o tempo de interação zero, o que torna impossível a descrição do problema a partir de uma integração de movimento simples. O método utilizado, a ser explicitado aqui, que resolve este problema é o evento dirigido.

Evento dirigido

A ideia do método para resolver o problema do força infinita é, ao invés de avançar o sistema em pequenos passos de tempo , avançar a simulação conforme as colisões forem ocorrendo. Para isso deve-se encontrar o par de partículas que colidirá no menor intervalo de tempo entre todas as partículas, denotaremos tal intervalo por , e, então, avançar o sistema. Neste ponto teremos dois objetos colados, portanto aqui deve ser feita a mudança de velocidades de tal forma a respeitar uma colisão elástica.

Fluxograma de um programa simples usando o método de evento dirigido com a otimização aqui explicitada.

Determinação do tempo de colisão

Os objetos a serem usados para o cálculo do tempo de colisão entre um par de partículas serão discos de raio , , cuja soma denotamos por . Portanto, segue que a condição de colisão é:

Com sendo o vetor posição da partícula e o tempo de colisão entre as partículas . Tal condição nos leva a determinação de a partir da expressão:

Onde , e .
Com isso, consegue-se determinar o valor de encontrando o menor valor de .

Mudança de velocidade em uma colisão elástica

Para fazer a mudança de velocidades temos que considerar o caso de colisão elástica entre as partículas , sendo impulso dado por:

.

Assim, a variação de velocidades pode ser determinada por:

e .

Otimização básica

Dado uma simulação de partículas, determinar o valor de é uma operação de ordem , ou seja, evitar fazer esse processo todo passo de tempo economiza grande parte do tempo computacional. Uma forma simples de fazer isso é determinar e armazenar o menor para cada partícula e o índice antes do loop temporal e em todo passo de tempo determinar o valor de a partir dos armazenados. Assim, a cada passo de tempo seria necessário apenas atualizar o valor de das partículas e das que colidiriam com ou .

Implementação computacional

Segue a implementação computacional, na linguagem de programação C, a função utilizada para o cálculo de . Importante notar que na formulação da funções-exemplo foram utilizadas condições de contorno periódicas no sistema.

void calc_dt_ij(double *xx, double *yy, double *vx, double *vy, double *dt_ij, int *colide, int i, int j){

//dt_ij é um vetor que armazena os menores valor de dt_ij para cada i
//colide é um vetor que armazena o valor de j correspondente a partícula que i colide no menor tempo dt_ij

  double delta_t_ij = INF;
  double dx, dy, dvx, dvy;
  double drdr, dvdv, drdv;
  double d;
  double X, Y;
  double sigma = 2*raio;
	
  X = xx[i] - xx[j];
  dx = fmod(X, Lx) - rint(fmod(X, Lx)/Lx)*Lx;
  Y = yy[i] - yy[j];
  dy = fmod(Y, Ly) - rint(fmod(Y, Ly)/Ly)*Ly;

  dvx = vx[i] - vx[j];
  dvy = vy[i] - vy[j];
  
  drdv = dx*dvx + dy*dvy;
  drdr = dx*dx + dy*dy;
  dvdv = dvx*dvx + dvy*dvy;
  
  d = pow(drdv, 2) - dvdv*(drdr - pow(sigma, 2));

  if(d > 0 && drdv < 0){
    delta_t_ij = (-1)*(drdv + sqrt(d))/(dvdv);
      if(delta_t_ij < dt_ij[i]){
	dt_ij[i] = delta_t_ij;
  	colide[i] = j;
      }
  }

}

Onde foram definidos INF como um número computacionalmente grande e raio como o raio das partículas, em particular consideramos os raios de todas como iguais. Para realizar a troca de velocidades do par de partículas que colidem de forma a ser coerente com uma colisão elástica, confeccionou-se a função abaixo.

void switch_veloc(double *xx, double *yy, double *vx, double *vy, int i, int j, double *Jx, double *Jy){

  double dx, dy, dvx, dvy;
  double drdr, dvdv, drdv;
  double deltavx, deltavy;
  double sigma = 2*raio, sigma2 = pow(sigma,2), d;
  double X, Y;
	
  X = xx[i] - xx[j];
  dx = fmod(X, Lx) - rint(fmod(X, Lx)/Lx)*Lx;
  Y = yy[i] - yy[j];
  dy = fmod(Y, Ly) - rint(fmod(Y, Ly)/Ly)*Ly;

  dvx = vx[i] - vx[j];
  dvy = vy[i] - vy[j];
  
  drdv = dx*dvx + dy*dvy;
  drdr = dx*dx + dy*dy;
  dvdv = dvx*dvx + dvy*dvy;

  Jx[i] = drdv*dx/sigma2;
  Jy[i] = drdv*dy/sigma2;

  deltavx = Jx[i];
  deltavy = Jy[i];

  //TROCA VELOCIDADES DA PARTICULA I
  vx[i] -= deltavx;
  vy[i] -= deltavy;

  //TROCA VELOCIDADES DA PARTICULA J
  vx[j] += deltavx;
  vy[j] += deltavy;

}

Após alterar as velocidades basta avançar o sistema em . Segue um exemplo de função.

void pos(double *xx, double *yy, double *vx, double *vy, double dt_min){

  int i;

  for(i = 0; i < (NP); i++){
    xx[i] = xx[i] + vx[i]*dt_min;
    yy[i] = yy[i] + vy[i]*dt_min;
  }

}
Plot exemplo de resultado de simulação. O par de partículas destacadas em vermelho são as colisoras.

Seguindo o fluxograma apresentado e utilizando as funções disponíveis chega-se, por exemplo, na animação ao lado, que mostra o decorrer da simulação conforme ocorrem as colisões.
















Adição do campo gravitacional

A adição de um campo gravitacional uniforme nesta simulação é bastante simples, basta definir uma nova variável, a qual aqui será chamada de g, e atribuir um valor de sua escolha a ela.
Vale lembrar que campo gravitacional aqui computado gera uma força central, ou seja, é radial em relação a cada partícula do sistema, de modo que só há força atuando na direção , no sentido negativo de .
Também é importante lembrar que aqui é desconsiderado a interação gravitacional entre quaisquer duas partículas, uma vez que assume-se que massas são pequenas e desprezíveis em frente a massa que forma o campo em que estão imersas.

Cálculo do tempo de colisão

Uma preocupação que alguém pode ter com a adição de um campo gravitacional na simulação seria de que seria necessário mudar o cálculo do tempo mínimo entre as colisões. No entanto, como será mostrado a seguir, isso não é necessário e o cálculo permanece o mesmo.
De maneira geral, sem cálculos, o que deve ser feito é abrir cada vetor de deslocamento, adicionando um termo de tempo quadrádico ligado à aceleração gravitacional do campo escolhido, e elevar ao quadrado a condição de que ocorra a colisão. A partir daí tudo que se é necessário fazer é trabalhar a álgebra, que no fim, eliminará todos os termos com a aceleração. Para verificar o cálculo completo veja na página Cálculo do tempo de colisão com aceleração.

Implementação

Como este campo só atua na componente , basta modificar as linha de código da seção anterior para adicionar um termo de tempo quadrático à variação da posição assim a linha de código fica:

yy[i] = yy[i] + vy[i]*delta_t - g*delta_t*delta_t/2;

E também adiciona-se uma aceleração em nesta mesma parte do programa, de modo que haverá uma atualização na velocidade em de , ficando:

vy[i] -= g*delta_t

Mais abaixo é possível ver uma animação de como fica o sistema com essas mudanças feitas.

Plot exemplo de resultado de simulação com um campo gravitacional uniforme de g=2. O par de partículas destacadas em vermelho são as que colidiram.