Estudo de reações químicas homogêneas via método de Monte Carlo

Study of homogeneous chemical reactions using Monte Carlo simulation

Resumo

By using the Monte Carlo simulation platform with probabilistic mathematical functions of the Boltzmann type, <img src="/img/revistas/qn/v36n5/a20res.jpg" align="absmiddle"/>, having activation energy and temperature as parameters, it was possible to assess important dynamic aspects of homogeneous chemical reactions of the types A → B and A <img src="/img/revistas/qn/v36n5/20s01.jpg" align="absmiddle"/> B. The protocol proved a useful tool in work with the basic concepts of Kinetics and Thermodynamics allowing its application both in class activities and for assisting experimental procedures.

Monte Carlo simulation; chemical reactions; chemical equilibrium

Monte Carlo simulation; chemical reactions; chemical equilibrium

EDUCAÇÃO

Estudo de reações químicas homogêneas via método de Monte Carlo

Study of homogeneous chemical reactions using Monte Carlo simulation

Ravir Rodrigues Farias; Luiz Augusto Martins Cardoso^*; Nemesio Matos Oliveira-Neto; Baraquizio Braga Nascimento Junior

Departamento de Química e Exatas, Universidade Estadual do Sudoeste da Bahia, Campus de Jequié, Rua José Moreira Sobrinho, s/n, 45206-190 Jequié - Bahia, Brasil

RESUMO

By using the Monte Carlo simulation platform with probabilistic mathematical functions of the Boltzmann type, , having activation energy and temperature as parameters, it was possible to assess important dynamic aspects of homogeneous chemical reactions of the types A → B and A B. The protocol proved a useful tool in work with the basic concepts of Kinetics and Thermodynamics allowing its application both in class activities and for assisting experimental procedures.

Keywords: Monte Carlo simulation; chemical reactions; chemical equilibrium.

INTRODUÇÃO

Simulação computacional, ou "experimento computacional", desempenha um papel fundamental na Ciência moderna, sendo considerada o terceiro ramo complementar às abordagens teóricas e experimentais tradicionais.^1,2 Tais experimentos computacionais têm sido amplamente aplicados como estratégias didáticas no ensino de Ciências (Física, Química, Biologia).^3-9 Uma das vantagens em se utilizar uma abordagem computacional pode estar no fato de que nem sempre conseguimos no laboratório, situações e/ou regimes extrapolados, tais como altas temperaturas e pressões. Também, do ponto de vista teórico, nem sempre podemos resolver analiticamente as equações que descrevem um sistema químico em estudo, dada sua complexidade; tais obstáculos podem ser facilmente contornados utilizando-se um computador. Outro fato que revela uma vantagem do método estocástico aqui apresentado é que outras abordagens, tais como a solução de equações diferenciais, necessitam de um conhecimento matemático que não é estudado em nível básico nem nos primeiros semestres de um curso de graduação em Física, Química e Biologia.

Além disso, de maneira geral, o uso de novas tecnologias, associado muitas vezes à experimentação, tem possibilitado a elaboração de aulas mais frutíferas em sala sendo utilizadas como ferramentas didáticas para controle de experimentos, na coleta e análise de dados experimentais, na simulação de fenômenos (físicos, químicos, biológicos) e na instrução dirigida.¹⁰ Em particular, nas simulações de fenômenos, os estudantes podem repetir diversas vezes o "experimento", tendo assim um maior tempo para entender, compreender e interpretar os resultados, além de verificar limites de validade dos modelos em questão.⁴

Nesse âmbito, propomos um modelo matemático para simular reações químicas homogêneas elementares do tipo A → B e AD B baseado no Método de Monte Carlo (MMC) o qual leva em consideração a temperatura em que a reação está ocorrendo, bem como as energias de ativação direta e inversa do sistema como parâmetros-chave para as simulações. Com esse modelo conseguimos reproduzir resultados importantes da Cinética Química e da Termodinâmica, os quais são exibidos e discutidos nas seções que seguem. Tal modelo pode ser usado como estratégia metodológica em todos os níveis de ensino, além da aplicação em pesquisa científica no ramo da Físico-Química. Por fim, ressaltamos que essa é a primeira etapa para a confecção de um software voltado para o ensino e a pesquisa de reações químicas. Já existem trabalhos que utilizam o MMC para o estudo cinético de reações químicas,^3,8 entretanto não consideram explicitamente as energias de ativação e a temperatura como parâmetros nas simulações, os quais julgamos essenciais para as discussões que seguem.

MODELO MATEMÁTICO

Neste modelo matemático consideramos um sistema que reage quimicamente, formado por N moléculas dos tipos A e B, de tal maneira que N = NA+ NB, sendo NA o número de moléculas do tipo A e NB o número de moléculas do tipo B. A fração em mols de cada espécie é definida como XA= NA/N e XB= NB/N. A cada instante é permitido às moléculas reagirem, isto é, A se transformar em B (reação direta) e B se transformar em A (reação inversa). Assumimos que as reações direta e inversa ocorrem com as seguintes probabilidades: ²

sendo kB a constante de Boltzmann, T a temperatura, EA e EB as energias de ativação das reações direta e inversa, respectivamente. Assim sendo, a dinâmica é baseada em parâmetros físico-químicos chave para a análise das reações químicas, a saber: temperatura e energias de ativação.

É de se esperar que as frações molares dos reagentes e produtos convirjam para aqueles valores previstos pela Termodinâmica, uma vez que uma das partes da dinâmica está baseada na atualização de Metropolis^1,2 dada pela Equação 1, que define as probabilidades de transição entre estados acessíveis ao sistema, fazendo necessariamente que o sistema relaxe para o equilíbrio, cujas distribuições de probabilidades tendem para distribuições exponenciais do tipo Boltzmann, o que, por exemplo, nos permite obter a Lei de Arrenhius.

Por simplicidade, definimos os parâmetros adimensionais q ≡ kBT/EA e EBA ≡ EB/EA . O primeiro é um parâmetro adimensional, relacionando a energia térmica kBT com a energia de ativação direta EA; o outro, a relação da energia de ativação inversa EB medida em relação à energia de ativação direta. Todas as análises serão feitas variando-se esses parâmetros.

Note que quando EB → ∞, PBA → 0, implicando numa reação irreversível. Assim sendo, para esse tipo de reação a dinâmcia é definida apenas pela Equação 1a.

Para realizar as simulações inicialmente indexamos com i = 1, 2, 3... N, todas as moléculas,³ depois seguimos os seguintes passos: (I) tentamos selecionar uma molécula Ai sorteando um número aleatório e uniformemente distribuído em [1,N]. Se não conseguimos selecionar uma molécula tipo Ai nessa tentativa, passamos para o passo (II). Se conseguirmos, um número aleatório z uniformemente distribuído em [0,1] é gerado e comparado com PAB. Se z < PAB, Ai se transforma em Bi. Caso contrário, Ai permanece como está. (II) Tentamos selecionar uma molécula do tipo Bi da mesma forma que na escolha de Ai. Não conseguindo, o passo finaliza. Se conseguimos, um número aleatório z uniformemente distribuído em [0,1] é gerado e comparado com PBA. Se z < PAB, Bi se transforma em Ai. Caso contrário Bi permanence como está.

O sorteio definido nos passos I e II, que nada mais é do que a urna de Ehrenfest, traduz as equações diferenciais envolvidas e dá a evolução da solução em função do tempo e, portanto, o comportamento das equações cinéticas.

Executando os passos I e II N vezes, sendo N o número total de moléculas, definimos o passo de Monte Carlo (MC) das simulações, de forma que em cada passo de MC atualizamos N moléculas. Para obter a evolução das frações em mols em função do passo de MC, preparamos nosso sistema inicialmente possuindo somente moléculas do tipo A, isto é, NA= N e NB= 0. Nos resultados exibidos aqui utilizamos N = 10⁴. A partir dessa condição inicial, realizamos as simulações e analisamos o comportamento de XA e XB em função do passo de MC (equivalente ao tempo), para valores pré-estabelecidos de q e EBA. Ao final realizamos uma média sobre 10² amostras. A necessidade de se realizar uma média sobre várias amostras se deve ao fato de que nos passos I e II, utiliza-se uma sequência de números "aleatórios". De fato, essa sequência de números é pseudoaleatória, uma vez que nenhum gerador de números produz verdadeiramente uma sequência de números aleatórios. Assim, realizamos uma média sobre 10² amostras para diminuir ao máximo a dependência dos resultados com a sequência de números gerada.

Utilizamos essas quantidades de molécula, 10⁴, e essa amostragem 10², pois com esses valores conseguimos resultados que se aproximam muito bem daqueles previstos pela Termodinâmica. A dependência dos resultados simulados com N e com o número de amostras será discutida mais à frente.

Salientamos que para os processos não reversíveis aqui estudados, excluímos o passo (II) pois, neste caso, uma molécula do tipo B nunca se transforma em uma do tipo IA. Um ponto que merece destaque é que a dinâmica definida pela Equação 1 e pelos passos I e II, que definem o passo de Monte Carlo, considera que a orientação relativa das moléculas nas colisões é sempre favorável à reação química, implicando em um fator estérico p = 1. Embora esse fator não esteja explicitamente incluído na dinâmica, nota-se que é possível obter um resultado muito forte que é a dependência da constante de velocidade com a temperatura e com a energia de ativação, isto é, a Lei de Arrenhius, que é nosso objetivo principal. Uma análise minuciosa acerca do fator estérico será realizada e esperamos em breve reportá-la.

Aqui ressaltamos que o programa foi escrito na linguagem de programação FORTRAN e compilado com o Intel^® FORTRAN Compiler edição 11.1 baixado de www.intel.com na versão não comercial. O programa não está ainda disponível e esperamos que em breve consigamos produzir um software livre e gratuito a partir deste ensaio.

RESULTADOS E DISCUSSÃO

Inicialmente realizamos experimentos computacionais para a reação do tipo A → B. Na Figura 1 mostramos XA em função do passo de MC, para quatro valores diferentes de q, isto é, diferentes valores de kBT/EA. Observa-se que XA diminui com o tempo até que todo reagente A seja consumido, restando apenas o produto (moléculas do tipo B). Nota-se que XA diminui mais rapidamente com o aumento de q. Podemos entender essa velocidade maior no consumo de A de duas maneiras distintas, porém equivalentes: uma, é imaginarmos que EA seja fixo e kBT esteja aumentando. Nesse caso o aumento da temperatura faz com que as colisões moleculares sejam mais efetivas, uma vez que mais moléculas do tipo A possuem energia cinética média suficiente para superar a energia de ativiação EA e se transformar em B. Assim sendo, a Figura 1 representa o consumo do reagente A em diferentes temperaturas. Outra possibilidade é pensarmos que EA esteja diminuindo e kBT seja fixo (o que faz q aumentar). Nesse caso, também observaríamos um consumo mais rápido do reagente A em função do tempo, pois quanto menor a energia de ativação, maior a chance de uma molécula do tipo A superar tal barreira energética. Nesse caso, poderíamos pensar em sistemas diferentes, pois temos diferentes valores de EA, mas todos à mesma temperatura. Um resultado importante e interessante que nosso modelo produz é que diferentes sistemas (diferentes EA) em temperaturas diferentes, isto é, kBT, possuem mesmo comportamento qualitativo e quantitativo bastando apenas possuírem mesmos valores de q = kBT/EA .

Outra característicca importante é mostrada no gráfico interno da Figura 1, onde se observa ln XA em função do passo de MC. Observa-se que XA decai exponencialmente e que o coeficiente angular aumenta com o aumento da temperatura. Do ajuste nas simulações obtivemos a expressão:

Os valores dos parâmetros X0A e k para cada valor de q, e o valor do erro padrão de k são exibidos na Tabela 1. Para o cálculo do erro padrão foi utilizado o software Grace-5.1.22 baixado livremente na plataforma linux para a distribuição Debian.

Como podemos observar, as simulações reproduzem bem os resultados reportados na literatura,^11,12 bem como o valor esperado de X0A póximo de 1, uma vez que iniciamos nossas simulações com esse valor.

Esses resultados mostram o comportamento de uma reação de primeira ordem. Neste ponto é importante analisarmos as flutuações e convergência de XA em função do número de moléculas. Para esse fim, utilizamos q = 0,35 e dois valores diferentes do número de moléculas, isto é, N=10² e 10⁴. Para outros valores de q o comportamento é semelhante. Observamos (Figura 1S, material suplementar) que as flutuações diminuem à medida que N cresce, implicando que no limite termodinâmico essas flutuações não seriam detectáveis. Outro fato importante é que o tempo de relaxação, isto é, o tempo em que a fração molar vai a 0, é da ordem de 10² passos de Monte Carlo (MC) para as duas curvas. Relembrando que definimos o passo de MC como sendo a atualização de N moléculas e explicitando o fator multiplicativo por N, tais tempos assumem os valores 100*10² = 10⁴ passos e 100*10⁴ = 10⁶ passos para as simulações com N=10² e N=10⁴, respectivamente, sendo agora cada passo a atualização de uma molécula.

Agora analisaremos a dependência das flutuações e convergências em termos da amostragem utilizada. Observando o comportamento de XA em função do tempo para N=10² e 10⁴ e duas amostragens diferentes, isto é, 1 e 10 (Figura 2S, material suplementar), notamos para N= 10² que as flutuações diminuem à medida em que aumentamos a amostragem utilizada (Figura 2Sa, material suplementar), implicando que o aumento da mesma tem o mesmo efeito do aumento do número de moléculas N. Para um número de moléculas maior, N=10⁴, notamos que este efeito se torna menos relevante uma vez que as flutuações já eram pequenas para esse valor de N utilizado (Figura 2Sb, material suplementar).

Neste ponto, discutiremos a dependência dos resultados obtidos para os parâmetros mostrados na Tabela 1 em função do número de moléculas N e da amostragem. Para esse fim realizamos simulações para q = 0,28 e uma única amostra. Para outros valores de q o comportamento é similar. Dos resultados simulados (Tabela 1S, material suplementar) observamos que à medida que aumentamos o número de moléculas N, o erro padrão associado a k diminui e o valor de k tende àquele previsto teoricamente 0.02811566. Um comportamento semelhante é observado quando fixamos N=10⁴ (Tabela 2S, material suplementar) e aumentamos a amostragem. Notamos também, em ambos os casos, que os valores de X0A tendem ao valor 1, o que era de se esperar pois iniciamos nossas simulações com esse valor.

Um comportamento similar, em função do número de moléculas N e da amostragem, é apresentado nas flutuações e convergências quando tratamos as reações reversíveis através do modelo aqui apresentado.

Dos resultados mostrados na Figura 1 podemos obter outras relações cinéticas das reações químicas de primeira ordem, a saber: a meia-vida e a Lei de Arrenhius. Para mostrar que o sistema segue a Lei de Arrenhius, que descreve como a constante de velocidade k depende da temperatura T e da energia de ativação EA, medimos os coeficientes angulares das curvas do gráfico interno da Figura 1 para cada valor de q já que cada coeficiente angular corresponde a um dado valor da constante de velocidade k. Foram utilizados além dos valores constantes na Figura 1, os valores de 0,18; 0,19; 0,20; 0,25; 0,30; 0,40 e 0,50 para q, os quais não aparecem na Tabela 1 e nem na Figura 1, para melhor clareza. A Figura 2 mostra os resultados de ln k, em função do recíproco de q. Os dados simulados se ajustam bem à linha contínua cuja função é dada por:

sendo b = 1,00 ± 0,01, que também está de acordo com os resultados da literatura.^11,12

Outra característica importante é a meia-vida, t½, em função de k. Para cada valor de k, medimos o tempo em que XA chega à metade do seu valor inicial. Os resultados são exibidos na Figura 3, na qual traçamos o gráfico de t½versus o recíproco de k. Os dados simulados se ajustam bem à equação:

sendo a= 1,00018 e a = 0,692597 , o que também está em excelente acordo com a literatura. A Equação 4 está representada na Figura 3 pela linha contínua. Para esse ajuste o coeficiente de correlação foi de 0,99999.

Para uma reação de primeira ordem, a meia-vida não depende da quantidade de matéria inicial do reagente, como ficou evidenciado na simulação. De fato, na Figura 4 exibimos o comportamento de XA em função do passo de Monte Carlo (MC) para q = 1,0 e 4 valores iniciais diferentes das frações em mols de A (X0A): 1,0; 0,9; 0,8 e 0,7.

Nesta figura observamos que as meia-vidas para as quatro situações são iguais, com valor de t½= 1,89 indicado pela linha vertical pontilhada no eixo do passo de Monte Carlo (MC).

A partir de agora, exibiremos os resultados para as simulações reversíveis do tipo A

B . Para esse caso, há mais um parâmetro para variar, E_BA, além do q. Na Figura 5 exibimos os resultados para os comportamentos de X_A (símbolos não preenchidos) e X_B (símbolos preenchidos) em função do passo de Monte Carlo (tempo). Na Figua 5a foi utilizado E_BA= 0,5 e q = 1,0 (diamantes) e 2 (quadrados). Na Figura 5b mostramos os resultados para os mesmos q's da Figua 5a e E_BA= 2,0. Observamos que para E_BA < 1 (E_B < E_A) o sistema possui, ao termalizar, uma fração em mols de reagente X_A maior do que a de B. Isso se deve ao fato da barreira de energia da reação inversa ser menor do que a da reação direta e, portanto, uma molécula do tipo IB possuir uma chance maior de superar tal barreira energética quando comparada com a chance de uma molécula do tipo A superar a barreira energética da reação direta. Uma análise análoga pode ser feita para o caso E_BA > 1 (E_B > E_A). O caso E_BA < 1 é chamado de reação "reagente favorecida" enquanto que E_BA > 1 reação "produto favorecida".

Outro aspecto importante resgatado em nossas simulações é a dependência das frações em mols finais com relação à temperatura. Em ambas as Figuras 5a e 5b exibimos as frações molares para duas temperaturas diferentes. Podemos notar na Figura 5a que o aumento da temperatura faz com que XA termalize em um valor menor e XB em um valor maior. Essencialmente esse resultado nos mostra dois aspectos importantes: um deles é a dependência da constante de equilíbrio, Keq = XB/XA, com a temperatura, uma vez que ambas as frações molares se alteram no equilíbrio. Outro aspecto é o princípio de Lê Chatelier aplicado ao efeito da temperatura. Observamos na Figura 5a que o aumento da temperatura desloca o sistema para a direção mais endotérmica. Uma análise similar pode ser feita na Figura 5b.

Na Figura 6, exibimos os resultados das frações molares em função do passo de Monte Carlo (tempo) para EBA = 1,0.

Neste caso observamos que as frações molares de equilíbrio são iguais a 0,5. O aumento da temperatura não as altera no equilíbrio, como é mostrado na Figura 6. Observa-se, nesse caso, que a constante de equilíbrio não depende da temperatura, devido ao fato das barreiras energéticas das reações direta e inversa serem iguais e, portanto, as moléculas do tipo A e B possuírem sempre a capacidade de superá-las com a mesma probabilidade e, como consequência, as quantidades das duas espécies de moléculas serão iguais para qualquer temperatura.

CONCLUSÕES E PERSPECTIVAS

Neste trabalho modelamos e simulamos reações químicas homogêneas elementares do tipo A → B e A

B. A modelagem foi realizada utilizando o Método de Monte Carlo (MMC) que inclui a temparatura e as energias de ativação como parâmetros fundamentais para o estudo em questão. Com este modelo simples pudemos recuperar todas as características básicas (ordem da reação, Lei de Arrenhius, princípio de Lê Chatelier) das reações químicas em estudo, principalmente aquelas relacionados com temperaturas e energias de ativação. Este modelo, portanto, apresenta grande potencial de aplicabilidade, também como aquele apresentado nas referências 3 e 8, podendo ser aplicado no estudo e na análise de sistemas químicos simples e complexos, auxiliando no âmbito didático-pedagógico em aulas de Físico-Química nos estudos de Termodinâmica e Cinética Química, possibilitando maior facilidade para a assimilação de conceitos dos processos moleculares subjacentes às transformações químicas em geral. Uma vantagem de simular tais sistemas em um computador se deve ao fato de que hoje em dia as escolas carecerem de laboratórios para o estudo experimental de reações químicas. Outra, é que o formalismo matemático para o estudo analítico via equações diferenciais é muito complicado, em todos os níveis, uma vez que a maioria dos estudantes não possui conhecimento suficiente para tal.

Por fim, ressaltamos que este trabalho pode ser muito útil para a futura confecção de softwares computacionais para o estudo de reações químicas mais complexas.

MATERIAL SUPLEMENTAR

Estão disponíveis em formato PDS, em http://quimicanova.sbq.org.br, os gráficos da evolução temporal de XA para uma amostragem e a comparação da evolução temporal em amostras diferentes. Constam também tabelas com dados do efeito do número de moléculas e do número de amostragens nos parâmetros de ajustes estatísticos.

AGRADECIMENTOS

Ao apoio parcial dado pela UESB, CNPq e CAPES e ao Prof. R. V. T. de Albuquerque pelas discussões proveitosas.

REFERÊNCIAS

1. Binder, K.; Heermann, D. W.; Monte Carlo Simulation and Statistical Physics, 2^nd ed., Springer-Verlag: New York, 1992.

2. Farias, R. R.; Cardoso, L. A. M.; Oliveira-Neto, N. M.; Exatas Online 2011, 2, 7.

3. López-Castillo, A.; Souza - Filho, J. C.; Quim. Nova 2007, 30, 1759.

4. Araújo, M. S. T.; Abib, M. L. V. S.; Rev. Bras. Ensino Fis. 2003, 25, 176.

5. Pereira, A. S. T.; Sampaio, F. F.; Ciência & Cognição 2008, 13, 51.

6. Eichler, M.; Del Pino, J. C.; Química. Nova na Escola 2000, 11, 10.

7. Novo, J. B. M.; Dias Júnior, L. C.; Quim. Nova 2011, 34, 707.

8. Winnischofer, H.; Araújo, M. P.; Dias Júnior, L. C.; Novo, J. B. M.; Quim. Nova 2010, 33, 225.

9. Hollauer, E.; Quim. Nova 1988, 11, 186.

10. Rosa, P. R. S.; Rev. Bras. Ensino Fis. 1995, 17, 182.

11. Panchenkov, G. M.; Lebedev, V. P.; Chemical Kinetics and Catalysis, 1^st ed., Mir Publishers: Moscow, 1976.

12. Atkins, P. W.; Physical Chemistry, 5^th ed., Oxford University Press: Oxford, 1994.

Recebido em 28/8/12

Aceito em 30/10/12

Publicado na web em 12/3/13

Datas de Publicação

• Publicação nesta coleção
  18 Jul 2013
• Data do Fascículo
  2013

Histórico

• Recebido
  28 Ago 2012
• Aceito
  30 Out 2012