ESTUDO CINÉTICO DE DECOMPOSIÇÃO TÉRMICA DE ESPUMAS RÍGIDAS DE POLIURETANO POR REDE NEURAL ARTIFICIAL

KINETIC AND THERMAL DECOMPOSITION STUDIES OF RIGID POLYURETHANE FOAMS MODELED BY ARTIFICIAL NEURAL NETWORKS

Resumo

Kinetic models of solid thermal decomposition are traditionally used for individual fit of isothermal experimental data. However, this methodology presents unacceptable errors in some regions of the data. To solve this problem, a neural network was adopted in this work. The implemented algorithm uses the rate constants as predetermined weights between the input and intermediate layer and kinetic models as activation functions of neurons in the hidden layer. The contribution of each model in the overall fit of experimental data is calculated as the weights between the intermediate and output layer. In this way, the phenomenon is better described as a sum of kinetic processes. Two rigid polyurethane foam samples: loaded with Al₂O₃ and no inorganic filler were used in this work. The R3 and D2 models described the thermal decomposition kinetic process for all temperatures for both foams with smaller residual error. However, the network, combining the kinetic models, presented residual errors on average 10² times lower compared to these individual models. The determined activation energy is 12.44 kJ mol^-1 higher for the loaded foam. This result corroborates the use of this material as flame retardant, even with the presence of a small amount of charge in its structure.

Keywords:
thermal decomposition; rigid polyurethane foams; polyurethane; artificial multilayer neural network

Keywords:
thermal decomposition; rigid polyurethane foams; polyurethane; artificial multilayer neural network

INTRODUÇÃO

Estudos cinéticos de decomposição térmica no estado sólido são de grande interesse científico e industrial.¹ Três metodologias relevantes no tratamento teórico de cinética de decomposição em sólidos são: (i) descrição do processo por modelos cinéticos,² (ii) associação e correção dos modelos cinéticos por rede neural artificial^3,4 e (iii) modelo isoconversional.⁵

A metodologia de ajuste por modelos cinéticos considera que as reações de decomposição térmica de sólidos ocorrem na interface do produto-reagente. Este processo cinético, com base na formação e crescimento de núcleos, pode ser estudado por análise de dados de decomposição térmica nos quais a redução da massa é medida num intervalo de tempo à temperatura constante. Os núcleos de reação são preferencialmente formados em imperfeições da estrutura e os modelos cinéticos são usados para explicar as isotermas experimentais, em que a fração de decomposição é medida ao longo do tempo. Um determinado modelo é escolhido devido a sua melhor correlação para ajustar os dados experimentais. No entanto, em muitas situações, o erro residual do modelo não é aceitável na descrição do processo total, embora, para regiões específicas de decomposição, o modelo possa ser apropriado, sendo apenas uma primeira aproximação. Os modelos podem ser representados por uma equação geral do tipo:²

(1) $\frac{\mathit{d}\left(\alpha \right)}{\mathit{dt}}=\mathit{k}{\left(\alpha \right)}^{1-\mathit{p}}{\left(1-\alpha \right)}^{1-\mathit{q}}$

em que α representa a fração de decomposição, t o tempo, k a constante de velocidade. Os parâmetros de ajuste p e q definem o modelo físico.

Em alguns trabalhos já publicados pelo nosso grupo, como pode ser verificado na literatura,^3,4,6 foi construído um programa computacional baseado na teoria de redes neurais de múltiplas camadas para determinar o conjunto de modelos que melhor descrevem a cinética de decomposição térmica de sólidos. Nestes trabalhos foram estudados a cinética de decomposição do acetato de Ródio e os fármacos Lamivudina e Efavirenz, usados no coquetel anti-HIV. Esta análise permite a correção matemática do modelo cinético e a determinação do processo físico global com erro experimental de ajuste das curvas bem menor que a utilização de apenas um dos modelos, o que consequentemente fornece informações relevantes sobre o material estudado. O método consiste na linearização da rede neural pela fixação dos pesos de interconexão entre a camada de entrada e a camada intermediária, fazendo com que cada neurônio da camada intermediária represente um modelo cinético a ser ponderado pelo neurônio na camada de saída. A rede neural representa um poderoso método de rotina para estudar reações de decomposição de sólidos.

A abordagem descrita por Ozawa consiste em tratar os materiais por técnicas termogravimétricas, em que a fração de decomposição é medida em função da temperatura em diferentes razões de aquecimento do sólido. Estas curvas experimentais são comparadas e os valores de decomposição, comuns a todas as curvas, são tratados a fim de se obter os parâmetros cinéticos (energia de ativação e fator de frequência) do processo. A equação geral para esta metodologia é:⁵

(2) $\ln \left(-\frac{\mathit{d}\left(1-\alpha \right)}{\mathit{dt}}\right)=-\frac{\mathit{E}}{\mathit{RT}}+\ln \left(\mathit{A}{\left(1-\alpha \right)}^{\mathit{n}}\right)$

No presente trabalho estaremos interessados em investigar a cinética de decomposição de espumas rígidas de poliuretano (ERPU), sem e com adição de carga inorgânica contendo 10% de Al₂O₃, pela metodologia de redes neurais. As espumas rígidas de poliuretano (ERPU) são uma das diversas formas em que se encontra o poliuretano (PU), um polímero formado essencialmente pela reação de condensação de um isocianato com um poliol. Nesta reação há dispersão de um gás, causada pelo aquecimento e evaporação de um agente de expansão durante o processo de polimerização, que origina uma estrutura tridimensional altamente reticulada de baixa densidade. As ERPU podem ser usadas em ampla faixa de temperatura, de -200 °C a +150 °C, por possuir alto ponto de amolecimento, além de ser resistente a produtos químicos e apresentar um baixo coeficiente de condutividade térmica, em virtude do gás de expansão preso nas unidades celulares e da baixa condutividade do próprio polímero. Devido a estas propriedades, é um dos mais eficientes isolantes térmicos usados na construção civil. Sua desvantagem é ser inflamável, mas esta característica pode ser contornada com a adição de substâncias que retardam a propagação das chamas, como dibromo propanol, o tris(dicloropropilfostato), (TCPP) e a alumina.^7-15

Em geral, os aditivos retardantes de chama são substâncias químicas que ao serem adicionadas aos polímeros interferirão nas condições da combustão. Essas substâncias podem ser classificadas de acordo com sua composição química como retardantes orgânicos reativos ou não reativos e os inorgânicos. Os dois primeiros apresentam compostos de fósforo e de halogênios, em especial cloro e bromo, e se diferenciam apenas pela presença ou não desses grupos funcionais presos quimicamente à cadeia polimérica. Dentre os retardantes de chamas inorgânicos podem ser citados o hidróxido de alumínio ou alumina trihidratata, material inerte de alta estabilidade. A alumina é preferível frente aos compostos contendo halogênio e fósforo que são extremamente tóxicos e prejudiciais ao meio ambiente.^16-19

O estudo cinético sobre a decomposição térmica das ERPU, sem e com adição de carga inorgânica, será realizado no presente trabalho com o principal objetivo de identificar a contribuição de modelos cinéticos na descrição dos dados experimentais de decomposição térmica e, consequentemente, identificar os mecanismos de decomposição das espumas, com e sem a presença de carga. O estudo sobre o mecanismo, nas duas situações, fornece informações valiosas sobre as propriedades fundamentais do material, como por exemplo suas propriedades mecânicas. Outra abordagem também a ser tratada é de investigar o acréscimo na energia de ativação referente ao evento de decomposição, provocado pela adição da carga inorgânica. Este acréscimo na energia é de fundamental importância para o uso a que se destina este material, pois terá sua atividade assegurada pela comprovada eficiência do retardante de chama a ele adicionado.

METODOLOGIA

Teoria dos modelos cinéticos

A cinética da decomposição térmica de materiais é baseada na nucleação e no crescimento dos núcleos ativos presentes na superfície do cristal.²⁰ A existência de pontos reativos separados, relacionados com as imperfeições dos cristais, provoca um aumento da energia de Gibbs do sistema e, consequentemente, a sua reatividade. Diversos modelos físicos são utilizados para descrever a fração decomposta, a, em função do tempo de decomposição. Sendo definido como a quantidade de perda de massa no tempo t normalizado para a massa total perdida.^20-22 Os modelos físicos estão apresentados na Tabela 1 e em geral são classificados de acordo com a forma da curva e com o tempo de aceleração e desaceleração.^2,20-24

Tabela 1
Modelos Cinéticos de Decomposição Térmica

Modelos	Símbolos	Equações Cinéticas
Aceleração
Lei Potencial	Pn	${\mathit{\alpha }}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathit{n}$}\right.}=\mathit{kt}+{\mathit{k}}_0,\mathit{n}=2,3,4\dots$
Sigmoide
Avrami-Erofeev	Am	${\left[-\mathit{ln}\left(1-\mathit{\alpha }\right)\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathit{m}$}\right.}=\mathit{kt}={\mathit{k}}_0,\mathit{m}=2,3,4\dots$
Avrami-Erofeev	Au	$\ln \frac{\alpha }{1-\alpha }=\mathit{kt}+{\mathit{k}}_0$
Prout-Tompkins	Ax	$\ln \frac{\alpha }{1-\alpha }=\mathit{k}\ln \mathit{t}+{\mathit{k}}_0,\mathit{k}>1$
Desaceleração
Modelo geométrico - contração
Contração Linear	R1	1-(1-α) = kt + k0
Contração Superficial	R2	$1-{\left(1-\alpha \right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}=\mathit{kt}+{\mathit{k}}_0$
Contração Volumétrica	R3	$1-{\left(1-\alpha \right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}=\mathit{kt}+{\mathit{k}}_0$
Modelo de Difusão
Uma Dimensão	D1	α2 = kt + k0
Duas Dimensões	D2	(1-α)ln(1-α) + α = kt + k0
Três Dimensões	D3	${\left[1-{\left(1-\alpha \right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\right]}^2=\mathit{kt}+{\mathit{k}}_0$
Ginstling-Brounshtein	D4	$1-\frac{2\alpha }{3}-{\left(1-\alpha \right)}^{\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}=\mathit{kt}+{\mathit{k}}_0$

As reações de decomposição, em que predomina a etapa de aceleração, caracterizam-se por apresentarem apenas o fenômeno de nucleação. Nessas reações a formação de núcleos pode ocorrer de forma instantânea ou com velocidade de nucleação constante e o modelo exponencial é mais usado para descrever o evento.^20-22 Já as reações em que se observa uma nucleação caótica seguida do crescimento desses núcleos, os modelos de Avrami-Erofeev ou Prout-Tompkins são mais utilizados.^23,24 Nas reações onde há apenas o crescimento dos núcleos, as curvas de desaceleração são mais apropriadas e dois fenômenos são responsáveis pela cinética dessas decomposições: a contração e a difusão.^20-24

A determinação do modelo físico no processo de decomposição de materiais é crucial para o estudo da cinética das reações e pode ser feita por análise de microscopia, através de softwares comerciais, como o Microcal Origin, mas que geralmente ficam restritos a no máximo três modelos para ajustar os dados experimentais ou por redes neurais artificiais, que ponderam a contribuição de todos os onze modelos físicos no processo.^3-4,25-28

Rede neural de multi camadas aplicada à decomposição térmica de sólidos

A metodologia adotada nesse trabalho utiliza um algoritmo proposto originalmente por Sebastião e colaboradores em 2003 e em 2004 para o estudo da decomposição térmica do acetato de ródio e acetado de ródio(II), respectivamente.^3,4 O software utilizado foi desenvolvido pelo grupo e está em processo inicial de registro. A rede multicamadas (MLP) aplicada a este tipo de problema possui três camadas, uma de entrada, uma intermediária e uma de saída. O diagrama da rede pode ser verificado na Figura 1.

Figura 1.
Arquitetura da rede MLP desenvolvida

Nas camadas de entrada e de saída há apenas um neurônio, representados pelos números 1 e 6 no diagrama, respectivamente. A camada de entrada, neurônio 1, recebe como entrada os dados experimentais de tempo, t. O neurônio 6, na camada de saída, gera como resposta da rede, Y_cal, dados de fração de decomposição, α(t). Na camada intermediária, o número de neurônios é variável, dependendo da quantidade de modelos cinéticos que se deseja incluir no estudo. No diagrama da Figura 1 são representados quatro neurônios nesta camada, 2, 3, 4 e 5. Os pesos de interconexão entre os neurônios seguem a nomenclatura, w_ij, em que i representa o neurônio da camada posterior e j o neurônio da camada anterior. Por exemplo, w₂₁ representa a interconexão entre os neurônios 2 e 1, sendo 2 o neurônio da camada posterior no sentido de propagação da informação na rede. Os pesos w_i0 indicam os bias, ou seja, são estímulos externos recebidos pelos neurônios, não provenientes de outros neurônios.

A metodologia adotada (Figura 1) por esta rede específica possui três etapas principais: a primeira etapa do algoritmo consiste em fixar os pesos de interconexão entre os neurônios da camada de entrada e a camada intermediária e fixar os bias dos neurônios na camada intermediária.^3,4,6 Essas constantes são armazenadas na matriz W₁, responsável pela linearização da rede. A matriz W₁ possui dimensão (n × 2), sendo n o número de modelos cinéticos considerados na rede. Na primeira coluna são alocados os pesos entre as camadas de entrada e intermediária, w_i1, e na segunda coluna são alocados os bias da camada intermediária, w_i0.

(3) ${\mathbf{W}}_1=\left(\begin{array}{cc}{\mathit{w}}_{21}& {\mathit{w}}_{20}\\ {\mathit{w}}_{31}& {\mathit{w}}_{30}\\ {\mathit{w}}_{41}& {\mathit{w}}_{40}\\ {\mathit{w}}_{51}& {\mathit{w}}_{50}\end{array}\right)$

Os valores w_i1 e w_i0 correspondem aos valores das constantes k e k₀ dos modelos cinéticos, nessa ordem, e são obtidos por um ajuste linear entre os modelos cinéticos e os dados experimentais de α(t) = Y_exp. No algoritmo geral das redes MLP, os pesos de interconexão entre todas as camadas são calculados após a minimização da função objetivo, descrita posteriormente pela Equação 8. As funções de ativação são funções suaves, crescentes e contínuas, e não possuem qualquer significado físico-químico relacionado ao problema específico. Na presente forma, os modelos cinéticos foram utilizados, após a fixação dos pesos da camada de entrada, como funções de ativação dos neurônios da camada intermediária, por apresentarem características similares. Esse algoritmo proposto descreve um diferente tipo de rede neural, em que a função ativação traz um significado físico-químico ao longo do processo de ajuste, ou seja, cada função ativação dos neurônios na camada intermediária é correspondente ao modelo cinético utilizado na rede para descrever o processo de decomposição térmica do sólido.

Os dados experimentais de tempo são recebidos pelo neurônio na camada de entrada como uma matrix do tipo:

$\mathbf{X}=\left(\begin{array}{ccccc}{\mathit{t}}_1& {\mathit{t}}_2& {\mathit{t}}_3& \dots & {\mathit{t}}_{\mathit{m}}\\ 1& 1& 1& \dots & 1\end{array}\right)$

Os neurônios na camada intermediária têm seus estados definidos após a multiplicação W₁X. Para a propagação da informação para a camada de saída, os neurônios precisam ser ativados. Então, na segunda etapa do algoritmo, é realizada uma transformação não-linear em que cada neurônio da camada intermediária é ativado pela função de ativação, f, correspondente ao modelo cinético. Como resultado tem-se a matriz B que determina os estados dos neurônios ativados, como mostra a Equação 4.^2-4

(4) $B=\mathit{f}\left(W_{1}X\right)$

Cada linha da matrix B corresponde aos dados de α(t) simulados por cada modelo cinético separadamente, como representados na Equação 5. Nesta etapa, α(t) é calculado analiticamente ou pela minimização do funcional que relaciona o modelo f(α) e as constantes k e k₀. Essa etapa é realizada por um ajuste não linear de mínimos quadrados (NLLS) com Levenberg-Marquardt.²⁹ Vale ressaltar que essa estrutura de rede não é rígida, sendo possível adicionar ou excluir modelos da camada intermediária, dependendo do sistema em estudo.^3-4,6

(5) $\mathbf{B}=\left(\begin{array}{c}{\alpha }_{{\mathit{D}}_1}\\ {\alpha }_{{\mathit{D}}_2}\\ {\alpha }_{{\mathit{D}}_3}\\ {\alpha }_{{\mathit{D}}_4}\\ {\alpha }_{{\mathit{R}}_1}\\ {\alpha }_{{\mathit{R}}_2}\\ {\alpha }_{{\mathit{R}}_3}\\ {\alpha }_{{\mathit{Am}}_4}\\ {\alpha }_{{\mathit{Am}}_2}\\ {\alpha }_{\mathit{AU}}\\ {\alpha }_{{\mathit{F}}_1}\end{array}\right)$

Os neurônios da camada intermediária são ponderados pelos pesos de interconexão entre a camada intermediária e a de saída pelo vetor W₂. Para o diagrama apresentado na Figura 1, o vetor W₂ pode ser exemplificado como apresenta a Equação 6.

(6) ${\mathbf{W}}_2=\left({\mathit{w}}_{62}{\mathit{w}}_{63}{\mathit{w}}_{64}{\mathit{w}}_{65}\right)$

E na terceira e última etapa ocorre uma transformação linear, resultando no vetor Y_cal que é a resposta da rede,

(7) $Y_{\mathrm{cal}}=W_2\mathit{f}\left(W_{1}X\right)$

Assim, a função erro da rede pode ser calculada pela diferença entre o estado do neurônio da camada de saída, e os dados de decomposição obtidos experimentalmente, Y_exp = α(t) como:

(8) $E={\left|Y_{\mathrm{cal}}-Y_{\exp }\right|}_2^2={\left|{\mathrm{BW}}_2-Y_{\exp }\right|}_2^2$

Com essa estrutura de rede é possível calcular a contribuição dos diversos modelos cinéticos através do algoritmo de otimização proposto por Levenberg-Marquart.²⁹ Neste algoritmo assume-se que em um mínimo, a derivada primeira da função deve ser zero. Para isto, a condição

$\frac{\partial \mathit{E}\left({\mathbf{W}}_2\right)}{\partial {\mathbf{W}}_2}=0$

faz com que a função erro seja minimizada e a expressão obtida a partir da Equação (8) após a aplicação desta restrição é:

$B^T{\mathrm{BW}}_2=B^T{Y}_{\exp }$

Desta forma, W₂ é calculado por:

(9) $W_2={\left(B^{T}B\right)}^{-1}{B}^T{Y}_{\mathit{exp}}$

A obtenção dos pesos de interconexão entre a camada de saída e a camada intermediária pelo vetor W₂ é uma metodologia particular da arquitetura de rede neural proposta neste trabalho, que foi estabelecida especificamente para este tipo de problema. É uma metodologia completamente diferente das redes neurais usuais, que utilizam, por exemplo, o algoritmo backpropagation para a obtenção dos pesos de interconexão da rede.

Regularização da solução da rede

Problemas que envolvem matrizes mal condicionadas ou com elevado valor de posto não são resolvidos através da Equação (9), pois as soluções não são satisfatórias, ou até mesmo incorretas. Para contornar este problema, a solução procurada deve ser aquela que apresenta a norma mínima do resíduo e também ser a solução de norma mínima. As soluções destes problemas mal condicionados devem ser obtidas pela remoção da singularidade através de um procedimento de regularização. No presente trabalho a regularização de Tikhonov foi utilizada.^29,30 Para isto, minimiza-se o funcional

(10) $E={\left|{\mathrm{BW}}_2-Y_{\exp }\right|}_2^2+\mathit{\lambda }{\left|{\mathrm{TW}}_2\right|}_2^2$

em que é um parâmetro a ser determinado e um operador do tipo

(11) $T{W}_2=a_o{\left|W_2-\stackrel{─}{W_2}\right|}_2^2+a_1{H}_1+a_2{H}_2$

com W₂ sendo uma solução de referência, por exemplo a solução não regularizada, H₁ o operador de derivada primeira, comumente referido por Jacobiano, e H₂ o operador de derivada segunda, expresso na literatura como Hessiana. As constantes a₀, a₁ e a₂ são nulas ou unitárias, dependendo do grau de regularização desejado.

Aplicando a condição de minimização, $\frac{\partial \mathit{E}\left({\mathbf{W}}_2\right)}{\partial {\mathbf{W}}_2}=0$, ao funcional da equação (10), obtém-se a expressão

$\left(B^{T}B+\mathit{\lambda }\left(a_{o}I+a_1{H}_1+a_2{H}_2\right)\right)W_2=B^T{Y}_{\exp }+\mathit{\lambda }\stackrel{─}{W_2}$

com I sendo a matriz identidade, H₁ e H₂, expressos comumente na literatura^29,30 como

$\begin{array}{c}{\mathbf{H}}_1=\left(\begin{array}{ccccc}1& -1& & & \\ -1& 2& -1& & \\ & \ddots & \ddots & \ddots & \\ & & -1& 2& -1\\ & & & -1& 1\end{array}\right)\hfill \\ {\mathbf{H}}_2=\left(\begin{array}{ccccccc}1& -2& 1& & & & \\ -2& 5& -4& 1& & & \\ 1& -4& 6& -4& 1& & \\ & \ddots & \ddots & \ddots & \ddots & \ddots & \\ & & 1& -4& 6& -4& 1\\ & & & 1& -4& 5& -2\\ & & & & 1& -2& 1\end{array}\right)\hfill \end{array}$

$\frac{\partial \mathit{E}\left({\mathbf{W}}_2\right)}{\partial {\mathbf{W}}_2}=0$

A solução W₂ pode então ser obtida pela expressão³⁰

(12) $W_2={\left(B^{T}B+T^{T}T\right)}^{-1}{B}^T{Y}_{\exp }$

A minimização do funcional proporciona simultaneamente a minimização da norma do resíduo e da norma da solução. Essas duas contribuições para o funcional são ponderadas pelo parâmetro l, em uma curva denominada curva-L.^29,30

Este procedimento é realizado escolhendo-se um intervalo de λ em que o processo de minimização deverá ser executado. Para cada λ determina-se a norma do resíduo, ${\left|{\mathbf{BW}}_2-{\mathbf{Y}}_{\mathbf{exp}}\right|}_2^2$ e a norma da solução, ${\left|{\mathbf{W}}_2\right|}_2^2$, em que W₂ é determinado pela Equação (12), utilizando-se as matrizes Jacobiana e Hessiana. Os valores da norma do resíduo e da norma da solução são agrupados em um gráfico que apresenta a forma de L e a melhor solução é determinada pela concavidade da curva, pois neste ponto há um balanço entre os menores valores da norma do resíduo e da norma da solução. Desta forma, a concavidade da curva L determina a melhor solução e o λ mais apropriado.³⁰

A determinação da contribuição dos diversos modelos cinéticos no ajuste dos dados experimentais possibilita o estudo físico ao longo de todo o processo de decomposição. Por exemplo, pode-se avaliar a sequência de eventos e as etapas físicas do processo global através da contribuição mais significativa para o ajuste em cada intervalo da curva experimental. A determinação mais precisa dos modelos acarreta na obtenção mais acurada dos parâmetros cinéticos e, consequentemente, auxilia no procedimento de mitigação ou aceleração do fenômeno, conforme seja o interesse experimental.

PARTE EXPERIMENTAL

As amostras de espumas rígidas de poliuretano (ERPU) utilizadas nesse trabalho contém glicerol como agente reticulador. A primeira amostra é uma ERPU sem carga (ERPU-0%Al₂O₃) e a segunda amostra recebeu cerca de 10% de carga inorgânica contendo alumina (ERPU-10%Al₂O₃). A carga utilizada é proveniente do filtro eletrostático da etapa de calcinação do processo Bayer, fornecido por indústrias de produção de alumínio do Estado de Minas Gerais.³¹

Ambas as espumas foram sintetizadas utilizando-se poliol de origem vegetal, produzido através da modificação química (esterificação) de óleo extraído das sementes de Ricinus communis, popularmente conhecida no Brasil como Mamona.^32,33 Essas espumas foram submetidas à decomposição térmica com intuito de obter suas curvas termogravimétricas (TG) dinâmicas, mostradas na Figura 2. As curvas TG foram obtidas em uma termobalança TGA-50H SHIMADZU, em cadinho de alumina, atmosfera de ar sintético com fluxo de 100 mL min^-1 e razão de aquecimento de 10 °C/min.

Figura 2.
Curva TG dinâmica para amostras das espumas rígidas de poliuretano contendo glicerol (a) sem adição de carga e (b) com adição de 10% de rejeito industrial contendo alumina, respectivamente

Quatro curvas TG isotérmicas experimentais, referentes ao início do processo de degradação dos grupos uretano, para cada ERPU, foram obtidas nas mesmas condições das curvas TG dinâmicas e são apresentadas na Figura 3.

Figura 3.
Isotermas experimentais às temperaturas (*)220 °C, (○)228 °C, (□)230 °C e (▷)232 °C e o ajuste fornecido pela rede (-) para as amostras de espumas rígidas de poliuretano contendo glicerol como agente reticulador, (a) sem adição de carga e (b) com adição de 10% de alumina (ERPU-10%Al₂O₃)

As curvas isotérmicas experimentais, mostradas na Figura 3 foram usadas nesse trabalho para investigar a cinética de decomposição dos grupos uretanos presentes no polímero. O principal objetivo é identificar o mecanismo cinético de decomposição e quantificar a energia de ativação relacionada ao evento de decomposição nas duas amostras. A verificação de um acréscimo no valor da energia de ativação para a espuma com carga é uma informação valiosa que assegura sua utilização como isolante térmico devido á comprovada eficiência do retardante de chama a ele adicionado. Também serão investigados os mecanismos de decomposição das espumas, que permitem a verificação do processo físico experimental e por consequência a melhor forma de sua mitigação.

RESULTADOS E DISCUSSÕES

Curvas TG

As curvas TG para as duas amostras apresentadas na Figura 2 exibem perfis semelhantes, com dois estágios em seus processos de decomposição. A primeira etapa de perda de massa é relativa à degradação dos grupos uretanos que ocorre na faixa de 215 °C a 350 °C. Já o segundo estágio refere-se à degradação dos resíduos da primeira decomposição que ocorre na faixa de 400 °C a 620 °C para ambas as espumas. Não há diferença macroscópica significativa entre as curvas isotérmicas a baixas temperaturas, pois as temperaturas são escolhidas quando o grau de degradação ainda é baixo. À medida que a temperatura aumenta, verifica-se o maior efeito da carga que interfere no processo global de decomposição.^16,18

Assim, como o processo de decomposição dos grupos uretanos ocorre em temperatura mais baixa, em torno de 228 °C, essa foi escolhida por ser a primeira observação experimental, além de apresentar perda de massa significativa nessa temperatura, para as respectivas amostras.

A temperatura em que o processo de decomposição se inicia, para os grupos uretano, nas espumas sem e com adição de alumina é a mesma, uma vez que a adição de 10% de carga não altera a estabilidade térmica dos compósitos. Porém, testes de chama realizados nestas mesmas amostras mostram uma redução de cerca de 30% nos tempos de autoextinção da chama para as espumas com 10% de carga, quando comparadas às ERPU sem carga, o que evidencia o efeito retardante de chamas desta carga na matriz polimérica.^32,33

Ajuste dos dados pelos modelos cinéticos e pela rede neural MLP

Inicialmente as constantes de velocidade k e k₀ foram obtidas a partir de um ajuste linear de cada isoterma experimental pelos modelos físicos linearizados, apresentados na Tabela 1.^4,6,23-24 O erro residual de ajuste de cada isoterma pelos modelos individuais para as amostras de ERPU-0%Al₂O₃ e ERPU-10%Al₂O₃ são apresentados na Figura 4. A Figura deve ser entendida como uma representação dos modelos cinéticos, descritos na abcissa, que apresentam um valor de erro residual de ajuste para cada isoterma, representadas na Figura 4 como (*) 220 °C, (○) 228 °C, (□) 230 °C e (▷) 232 °C. Como pode ser observado na Figura 4a, os modelos D2 e R3 se mostraram mais adequados por apresentarem menores erros residuais na descrição das quatro isotermas estudadas. Na Figura 4b, para a amostra com 10% de carga, podemos notar que o modelo cinético D2 e o modelo de contração volumétrica, R3, são também os modelos cinéticos que melhor descrevem a decomposição térmica do grupo uretana presente na espuma em todas as temperaturas. Entretanto, para esta amostra, o modelo R2 também apresentou similar valor de erro de ajuste para os dados à temperatura de 220 ^oC.

Figura 4.
Erro residual de ajuste dos modelos cinéticos individuais, resultado do ajuste linear para as amostras de espumas rígidas de poliuretano, (a) sem adição de carga e (b) com adição de 10% de carga contendo alumina (ERPU-10%Al₂O₃) nas temperaturas (*)220°C, (○)228°C, (□)230° e (▷)232°C, respectivamente

Uma rede neural MLP foi proposta com o intuito de indicar uma função matemática que pondera a contribuição dos modelos físicos para melhor descrever o processo de decomposição térmica para as duas amostras nas quatro temperaturas. O estudo foi realizado utilizando-se uma rede neural com todos os modelos físicos sendo considerados como função de ativação na camada intermediária. Os dados de k e k₀ obtidos pelo ajuste dos modelos individuais aos dados experimentais foram organizados na matriz W₁ para serem utilizados no algoritmo da rede. A Figura 3 apresenta os ajustes promovidos pelas redes neurais às isotermas experimentais das espumas. Como pode ser observado, foi obtida uma excelente concordância entre os dados experimentais de decomposição térmica e o ajuste da rede neural. Os erros residuais obtidos pela rede são da ordem de 10^-5 para ambas as espumas analisadas, com a solução regularizada e sem regularização, conforme apresentado na Tabela 2.

Tabela 2
Norma do resíduo e norma da solução, sem regularização (Y_cal e com regularização (Y_{cal_reg}), para as espumas rígidas de poliuretano ERPU-0%Al₂O₃ e ERPU-10%Al₂O₃ em suas respectivas temperaturas de análise

	ERPU	Temperatura /ºC
		220	228	230	232
||Ycal - Yexp||22/10-5	0% Al2O3	1,538	3,492	6,069	2,127
	10% Al2O3	0,2929	0,1100	3,290	3,654
||Ycal||22	0% Al2O3 10% Al2O3	66,11 177,3	724,4 15,13	226,3 11,43	143,1 11,21
||Ycal_reg - Yexp||22/10-5	0% Al2O3 10% Al2O3	0,4595 0,6539	7,390 0,1304	0,1952 4,206	5,534 4,296
||Ycal_reg||22	0% Al2O3 10% Al2O3	0,3126 1,343	3,315 0,5350	2,668 0,2020	1,567 0,2360

A Figura 5 apresenta uma comparação do erro residual da rede MLP e os modelos D2 e R3. Os modelos D2 e R3 foram escolhidos devido ao melhor ajuste dos dados experimentais em todas as temperaturas em comparação com os outros modelos. Entretanto, a rede MLP reduziu substancialmente o erro residual de ajuste dos dados em todas as temperaturas para as duas amostras. Esta redução do erro residual da rede frente ao ajuste individual dos modelos pode ser explicada com base no fator de correção W₂, que introduz um grau de liberdade ao sistema e indica a contribuição de cada modelo no ajuste total. Em casos de modelos assintóticos, como os de Avrami-Erofeev e contração de volume, este parâmetro representa também o valor do resíduo remanescente do processo de decomposição. Por exemplo, ao isolar alfa na equação do modelo fornecido na Tabela 1 e considerando apenas o modelo de contração de volume, R3, na rede, o estado do neurônio na camada de saída pode ser expresso como y_calc = W₂ × [(kt + k₀ - 1)³ + 1], e no intervalo de tempo do experimento, o valor assintótico da curva de decomposição é corrigido pelo parâmetro W₂ adicionado pela rede, representando o resíduo final do material. Ou seja, se apenas o modelo R3 for considerado como neurônio na camada intermediária da rede neural, ainda assim este ajuste promovido pela rede será superior se comparado ao ajuste do modelo individual, pois um parâmetro físico foi introduzido ao modelo matemático.

Figura 5.
Erro residual das redes neurais MLP (□) e modelos D2(○) R3(*) para as amostras de espumas rígidas de poliuretano, (a) sem adição de carga contendo alumina (ERPU-0%Al₂O₃) e (b) com adição de 10% de carga contendo alumina (ERPU-10%Al₂O₃) nas temperaturas 220 °C, 228 °C, 230°C e 232 °C, respectivamente

Menores erros de ajustes dos dados experimentais representam maior exatidão no cálculo dos parâmetros cinéticos. A função matemática fornecida pela rede para a descrição dos dados experimentais ofereceu um ajuste aprimorado de no mínimo 10² vezes se comparado ao melhor modelo individual. Essa função determinada pela rede indica fundamentalmente que o processo global de decomposição não deve ser considerado como um processo de etapa única em toda sua extensão. Contrariamente, assume-se que cada região da curva deve ser descrita por um processo físico representado por um modelo cinético.

Uma outra abordagem possível é a utilização da função da rede por outros métodos de análise cinética, como por exemplo, o método de Ozawa, em que dados de temperatura devem ser determinados para o mesmo grau de conversão em curvas com diferentes taxas de aquecimento.⁵ Entretanto, para esta metodologia, deve-se considerar apenas um processo físico durante toda a decomposição, ou seja, $\raisebox{1ex}{$\mathit{d}\alpha $}\!\left/ \!\raisebox{-1ex}{$\mathit{dt}$}\right.=\mathit{k}\left(\mathit{T}\right)\mathit{f}\left(\alpha \right)$. Neste tratamento, a arquitetura da rede é definida por um neurônio na camada de entrada, um número variável de neurônios na camada intermediária e um neurônio na camada de saída e todos os pesos de interconexão entre os neurônios devem ser otimizados. Os dados de entrada são os dados de temperatura e a saída são os dados de fração de decomposição. Este tipo de rede fornece também bons resultados de ajuste, como pode ser verificado na literatura.^4,28

A Figura 6 apresenta os dados da contribuição de cada modelo, W₂, em cada temperatura para a rede com os onze modelos na camada intermediária. A solução apresentada é a solução regularizada, equação (12), com as constantes a₀, a₁ e a₂ consideradas como unitárias. A solução regularizada foi obtida a partir de um balanço da norma do resíduo e da norma da solução, ponderados pelo parâmetro λ. Como a matriz B na equação (12) é uma matriz mal condicionada, cond ≅ 10⁴, a regularização torna-se necessária para a obtenção da solução mais adequada. A Tabela 2 apresenta uma comparação entre a norma do resíduo e a norma da solução para a solução obtida sem regularização e a solução regularizada. A norma da solução regularizada é cerca de 100 a 200 vezes menor que a solução não regularizada e a norma do resíduo é da mesma ordem de grandeza, o que justifica a escolha da solução regularizada.

Figura 6.
Contribuição dos modelos cinéticos na função da rede neural para ajuste dos dados experimentais para a amostra sem carga (a) e com carga (b) nas temperaturas de análise (*) 220 oC, (o) 228 oC, (○)230 oC e (▷) 232 ^oC

Analisando a Figura 6 verifica-se um balanço entre os coeficientes da função da rede que proporcionam um melhor ajuste dos dados experimentais ponderando a norma do resíduo e da solução. Os valores dos coeficientes não relacionam a contribuição do modelo ao menor erro de ajuste individual dos dados experimentais para o processo global, mas informa a contribuição de cada um quando o processo é considerado em multietapas, sendo cada região do processo melhor descrita por um modelo.

Esta análise é interessante, pois a função matemática que descreve o estado do neurônio da camada de saída é composta por um somatório de onze funções correspondentes aos modelos cinéticos, ponderadas por sua contribuição no ajuste total da curva. Desta forma, nem sempre o modelo que apresenta um erro individual menor é o que tem maior contribuição quando os outros modelos são também considerados, pois cada modelo vai contribuir para determinadas regiões no ajuste total da curva.

Por exemplo, para a amostra sem carga à 220 ^oC, a Figura 6a apresenta o processo de decomposição pela rede neural como α(t) = 0,33D1(t) + 0,30D2(t) - 0,03D3(t) + 0,10D4(t) + 0,04R1(t) + 0,15R2(t) + 0,13R3(t) + 0,02Am4(t) - 0,04Am2(t) - 0,18Au(t) + 0,21F1(t). Esta função descreve a isoterma experimental com erro de ajuste 10³ vezes menor que o modelo D2, como pode ser verificado na Figura 5a. Este é o modelo individual que apresenta menor erro de ajuste nesta temperatura, utilizando o procedimento convencional. Esta função da rede pode ser então utilizada como um potente filtro ou interpolador/extrapolador para os dados experimentais, o que possibilita a determinação das constantes de velocidade com maior precisão. Após a determinação da função global da rede para cada isoterma, as constantes de velocidade são determinadas ajustando-se as funções dos modelos corrigidos aos dados experimentais. Este procedimento de correção do modelo cinético é realizado quando utiliza-se apenas um neurônio na camada intermediária da rede.

Considerando os modelos D2 e R3 como principais modelos para descrever o fenômeno de decomposição das espumas, assumimos que a decomposição se processa com o crescimento dos núcleos na superfície do material e os dois fenômenos de contração e difusão são responsáveis pela cinética da decomposição. A contração é responsável pelo desenvolvimento rápido de núcleo em toda a extensão da superfície do cristal. Já a difusão é responsável pelo controle da taxa de reação, uma vez que a continuidade da reação requer o transporte dos reagentes para a camada de produto.

A Figura 7 apresenta o ajuste da isoterma experimental à 220 ^oC pelos modelos D2 e R3, corrigidos pela rede. Ou seja, foi utilizada uma rede neural com apenas o modelo D2 e outra com o modelo R3, separadamente, como neurônio na camada intermediária, para o ajuste dos dados experimentais. Pode-se verificar as regiões onde os dois modelos melhor se ajustam aos dados experimentais. Na figura também é apresentado o ajuste promovido pela função da rede quando os onze modelos são considerados para descrever todo o processo. Neste ponto é possível destacar a superioridade desta função, que descreve o fenômeno em multietapas, frente aos modelos cinéticos individuais, mesmo sendo corrigidos pela rede.

Figura 7.
Isoterma experimental à 220 ^oC (o). A figura (a) apresenta o ajuste da função corrigida pelo modelo D2 (+) e a Figura (b) apresenta o ajuste pelo modelo corrigido R3(+). Em ambas as figuras o ajuste pela rede neural, considerando onze modelos cinéticos, é representado pela curva tracejada

Os parâmetros cinéticos de energia de ativação, Ea, e fator de frequência, A, foram calculados considerando-se a teoria de Arrhenius e as constantes de velocidade obtidas por todos os modelos cinéticos corrigidos, nas quatro temperaturas investigadas para as duas amostras. Os resultados são apresentados na Tabela 3. O gráfico de Arrhenius obtido para o modelo R3 é apresentado como exemplo na Figura 8, em que podemos observar que a linearidade entre os dados garante a validade da teoria. Todos os outros gráficos apresentaram comportamento similar.

Tabela 3
Energia de Ativação e o Fator de Frequência de acordo com os modelos cinéticos para as espumas ERPU-0%Al₂O₃ e ERPU-10%Al₂O₃

Modelo Cinético	ERPU-0%Al2O3		ERPU-10%Al2O3
	Ea /kJ mol-1	lnA	Ea /kJ mol-1	lnA
D1	52,27	3,296	63,64	6,154
D2	53,58	3,569	65,34	6,498
D3	52,12	2,668	65,54	5,960
D4	53,74	2,412	65,58	5,347
R1	49,97	2,533	59,49	4,984
R2	51,69	2,927	62,39	5,623
R3	51,63	2,773	63,12	5,643
Am4	49,82	2,456	62,05	5,540
Am2	49,67	3,138	62,50	6,340
AU	48,73	4,217	64,27	8,083
F1	48,32	3,709	64,50	7,675

Figura 8.
Grafico de Arrhenius (ln k x 1/T) para o modelo R3 para as amostras de espumas rígidas de poliuretano, (a) sem adição de carga e (b) com adição de 10% de carga respectivamente

Considerando os modelos R3 e D2 como os modelos que mais contribuíram para o menor erro residual da rede e são também os modelos de menores erros de ajuste individual, o valor da energia de ativação para a espuma ERPU-0%Al₂O₃ é de 51,63 kJ mol^-1 e 53,58 kJ mol^-1, nessa ordem. Estes valores são condizentes com a literatura,¹⁷ em que a espuma sem carga foi analisada e a energia de ativação determinada pelo método de Ozawa modificado como, aproximadamente, 50 kJ mol^-1 a baixas temperaturas.

Como pode ser verificado na Tabela 3, a energia de ativação do processo de decomposição para a espuma com carga é sempre maior que para a espuma sem carga, o que indica a necessidade de um pouco mais de energia no processo de decomposição dos grupos uretanos presente em sua estrutura, em média 12,44 kJ mol^-1. Esse fato corrobora com as propriedades retardantes de chama do material, mesmo com a presença de uma pequena quantidade de carga em sua estrutura. Essa carga afeta a cadeia polimérica da espuma por meio da formação de uma camada refratária que atua como uma barreira protetora.^16,18

CONCLUSÃO

O processo de decomposição térmica de espumas rígidas de poliuretano com glicerol como reticulador, com e sem adição de carga inorgânica, é estudado no presente trabalho. Uma rede neural artificial do tipo MLP foi proposta para o estudo, com a particularidade de assumir a função de ativação dos neurônios da camada intermediária como modelos cinéticos de decomposição.

Uma característica fundamental das redes neurais MLP é que os neurônios devem ser ativados por funções que tenham comportamento suavemente crescente e contínuo. As funções que descrevem os modelos físicos de decomposição de sólidos apresentam essa característica e por isso foram utilizados no presente trabalho como função ativação desta específica arquitetura de rede proposta. Por exemplo, em uma rede de três camadas os modelos podem ser utilizados se, e somente se, os pesos de interconexão entre a camada de entrada e a oculta forem fixados de acordo com as constantes cinéticas estabelecidas por cada modelo no ajuste dos dados experimentais. Este procedimento resulta na linearização da rede neural e faz com que as funções de ativação apresentem conceitos químicos. Em geral, as redes neurais MLP minimizam a função erro ajustando todos os pesos de interconexão entre os neurônios. Desta forma, nenhuma informação química é obtida por estas funções. Por outro lado, a linearização da rede por esta metodologia garante que informações químicas não serão perdidas no processo de ajuste pela rede neural, ou seja, a rede neural, sob esta concepção, não deve ser entendida como um simples mecanismo de ajuste matemático.

Esta rede MLP é uma poderosa ferramenta para tratar o fenômeno de decomposição térmica de sólidos, pois os erros residuais de ajuste são menores do que os resultados obtidos quando são utilizados os modelos individualmente. Além disso, com essa abordagem é possível calcular a contribuição de cada modelo na descrição experimental dos dados. Esta metodologia funciona, pois a rede neural propõe considerar o valor assintótico correto da fração de massa além de alterar a escala de tempo do processo nos modelos cinéticos. Essas alterações permitem a redução no erro residual em média 10² vezes, para ambas as espumas, quando comparado com o melhor modelo cinético individual que descreve o processo, R3 e D2.

Com este melhor ajuste, a rede neural MLP possibilita o estudo físico detalhado do processo e, consequentemente, o cálculo mais preciso dos parâmetros cinéticos como a energia de ativação e o fator de frequência. A metodologia apresentada nesse trabalho não é restrita às espumas rígidas estudadas e pode ser aplicado a outros sistemas, como inúmeros trabalhos abordados na literatura,^3,4,27-29 o que sugere um poderoso método de rotina no processo de decomposição térmica de sólidos.

AGRADECIMENTOS

Os autores agradecem as agências de fomento Fapeming, Capes, CNPq e PRPq/UFMG pelo apoio financeiro.

REFERÊNCIAS

• ¹
  Tseng, A.; Seet, J.; Phillips, E. J.; Brit. J. Clin. Pharmacol. 2014, 79, 182.
• ²
  Ng, W. L.; Austr. J. Chem. 1975, 28, 1169.
• ³
  Sebastião, R. C. O.; Braga, J. P.; Yoshida, M. I.; J. Therm. Anal. Calorim. 2003, 74, 811.
• ⁴
  Sebastião, R. C. O.; Braga, J. P.; Yoshida, M. I.; Thermochim. Acta 2004, 412, 107.
• ⁵
  Ozawa, T.; Bull. Chem. Soc. Jpn. 1965, 38, 1881.
• ⁶
  Ferreira, B. D. L.; Araujo, B. C. R.; Sebastião, R. C. O.; Yoshida, M. I.; Mussel, W. N.; Fialho, S. L.; Barbosa, J.; J. Therm. Anal. Calorim. 2017, 127, 577.
• ⁷
  Veronese, V. B.; Menger, R. K.; Forte, M. M. D. C.; Petzhold, C. L.; J. Appl. Polym. Sci 2011, 120, 530.
• ⁸
  Hentschel, T.; Münstedt, H.; Polymer 2011, 2, 3195.
• ⁹
  Demiryürek, O.; Koç, E.; Fibers Polym. 2009, 10, 237.
• ¹⁰
  Ashida, K.; Polyurethane and Related Foams: Chemistry and Technology, 6^th ed., CRC Taylor & Francis: Boca Raton, 2007.
• ¹¹
  Eaves, D.; Em Handbook of Polymer Foams, Eaves, D., ed.; Shawbury: Rapra Technology Limited, 2004, cap. 4.
• ¹²
  Chuayjuljit, S.; Maungchareon, A.; Savari, O.; J. Reinf. Plast. Compos. 2010, 29, 218.
• ¹³
  Oertel, G.; Polyurethane Handbook, 2^th ed., Hanser Gardner Publications: New York, 1994.
• ¹⁴
  Rabello, M.; Aditivação de Polímeros, 1° ed., Artliber Editora: São Paulo, 2000.
• ¹⁵
  Sato, T.; Thermochim. Acta 1985, 88, 69.
• ¹⁶
  Modesti, M.; Lorenxetti, A.; Simioni, F. and Camino, G.; Polym. Degrad. Stab. 2002, 77, 195.
• ¹⁷
  Li, K.; Pau, D. S.; Zhang, H.; Fire Mater. 2016, 40, 800.
• ¹⁸
  Oliveira, I. R.; Pandolfelli, V. C.; Cerâmica 2009, 55, 33.
• ¹⁹
  Hatakeyema, H.; Tanamachi, N.; Matsumura, H.; Hirose, S.; Hatakeyama, T.; Thermochim. Acta 2005, 431, 155.
• ²⁰
  Yeremin, F. N.; The Foundations of Chemical Kinetics, Mir: Moscow, 2001.
• ²¹
  Galwey, K.; Brown, M. E.; Stud. Phys. Theor. Chem. 1999, 86, 75.
• ²²
  Rodante, F.; Vecchio, S.; Tomassetti, M.; Thermochim. Acta 2002, 394, 7.
• ²³
  Tompkins, P. W. M.; Em Chemistry of the Solid State; Garner, W. E.; ed.; Butterworths Scientific Publication: London, 1983, p. 184.
• ²⁴
  Avrami, M.; J. Chem. Phys 1939, 7, 1103.
• ²⁵
  Tarrío-Saavedra, J.; López-Beceiro, J.; Naya, S.; Francisco-Fernández, M.; Artiaga, R. J.; J. Therm. Anal. Calorim. 2014, 118, 1253.
• ²⁶
  Straszko, J.; Biedunkiewicz, A.; Strzelczak, A.; Pol. J. Chem. Technol. 2008, 10, 21.
• ²⁷
  Isabel Moreno, A.; Font, R.; Conesa, J. A.; Fuel 2017, 192, 169.
• ²⁸
  Krawczyk, M.; Figiel, P.; J. Therm. Anal. Calorim. 2016, 126, 585.
• ²⁹
  Levenberg, K; Quart. Appl. Math. 1944, 2, 164; Marquardt, D. W.; SIAM J. Appl. Math. 1963, 11, 431.
• ³⁰
  Braga, J. P.; J. Math. Chem. 2001, 29, 151; Tikhonov, A. N.; Arsénine, V.; Méthods de Résolution de Problemes Mal Posés, Mir: Moscow, 1987.
• ³¹
  Yoshida, M. I.; Silva, V. R.; Pinto, P. C. C.; Sant'Anna, S. S.; Silva, M. C.; Carvalho, C. F.; J. Therm. Anal. Calorim. 2012, 109, 1429.
• ³²
  Silva, V. R.; Mosiewicki, M. A.; Yoshida, M. I.; Silva, M. C.; Stefani, P. M.; Marcovich, N. E.; Polym. Test. 2013, 32, 438.
• ³³
  Silva, V. R.; Mosiewicki, M. A.; Yoshida, M. I.; Silva, M. C.; Stefani, P. M.; Marcovich, N. E.; Polym. Test. 2013, 32, 665.

Datas de Publicação

• Publicação nesta coleção
  10 2017

Histórico

• Recebido
  31 Mar 2017
• Aceito
  21 Ago 2017