Avaliação de Diferentes Abordagens na Solução do Problema de Equilíbrio Sólido-Líquido em Óleos Parafínicos

RESUMO

Neste trabalho investigou-se a eficiência computacional na solução do problema de equilíbrio termodinâmico sólido-líquido (ESL) em óleos parafínicos. Foram avaliados diversos modelos termodinâmicos para a fase líquida: a equação de estado de Peng-Robinson (PR-EOS) e os modelos de coeficiente de atividade Wilson preditivo, UNIversal QUAsiChemical (UNIQUAC) preditivo e UNIQUAC Functionalgroup Activity Coefficient (UNIFAC), além do modelo de solução ideal. A fase sólida é caracterizada pelo modelo multisólido e a previsão de formação de fase sólida é inicialmente avaliada por um teste de estabilidade termodinâmica. Posteriormente, o sistema de equações não-lineares (SENL) formado pelas equações de equilíbrio termodinâmico e de balanço material é resolvido pelo método de Newton multivariável. Até onde pudemos verificar, não há estudos sistemáticos sobre o efeito do tipo de modelo (e sua complexidade) no tempo de computação necessário para a solução do problema. Por outro lado, tendo em conta o elevado número de componentes nas frações de petróleo estudadas, aspectos como a robustez (frente a diferentes estimativas iniciais) e tempo de computação são aspectos importantes neste contexto e que foram analisados neste trabalho.

Palavras-chave:
modelo multisólido; deposição de parafinas; equilíbrio sólido-líquido

ABSTRACT

In this paper, the computational efficiency in the solid-liquid equilibrium (SLE) problem for paraffinic oils was examined. Several different thermodynamic models were evaluated for liquid phase: the Peng-Robinson equation of state (PR-EOS), and the activity coefficient models for ideal solution, the predictive Wilson model, the predictive Universal Quasi-chemical (UNIQUAC) model, and UNIQUAC Functional-group Activity Coefficient (UNIFAC) model. The solid phase is characterized by the multisolid model and the prediction of solid phase formation is initially estimated by a thermodynamic stability test. Afterwards, the systems of non-linear equations (SNLE) formed by the thermodynamic equilibrium and the material balance equations was solved using Newton-Raphson method (multivariate). As far as we can see, there are no systematic studies addressing the effect of the model (and its complexity) on the computational effort of the algorithm. On the other hand, considering the high number of different components in the oil fractions studied, the robustness (due to different initial estimates), and the computational time are important aspects and were analyzed in this work.

Keywords:
multisolid model; paraffin deposition; solid-liquid equilibrium

1 INTRODUÇÃO

O problema do equilíbrio sólido-líquido tratado aqui consiste em: dada a temperatura, a pressão e a composição (global) de um petróleo parafínico, determinar a quantidade de cada componente na fase líquida e na fase sólida, bem como as quantidades das fases. A elucidação desta questão ocupa um papel chave na solução de problemas para estimar o crescimento da camada de parafina em um duto pelo qual escoa o petróleo. Durante o transporte deste óleo, do subsolo marítimo até a plataforma de petróleo, seus componentes mais pesados tendem a se solidificar e se aglutinar na superfície do duto que o transporta, diminuindo a seção transversal útil desse duto. Como resultado, há um aumento importante no custo de bombeamento do óleo, e, eventualmente, haverá a necessidade de se trocar o duto. Tendo em conta tal cenário, a simulação da deposição de parafinas vem sendo abordada na literatura especializada nos últimos anos. Para simular a precipitação de parafinas, comumente são utilizados modelos termodinâmicos que são classificados em:

• Modelo multisólido: modelo em que os precipitados são tratados como componentes puros e imiscíveis das demais fases sólidas ^{6), (10), (11}.

• Modelo de solução sólida: modelo em que os precipitados se dissolvem uns nos outros e formam uma única fase sólida ^{2), (3), (4), (18), (19}.

Geralmente, os problemas de escoamento de óleos parafínicos se baseiam na solução numérica de equações da Mecânica dos Fluidos e de Transferência de Calor e Massa, os quais permitem prever a quantidade de parafina depositada nos dutos. Para realizar a previsão da quantidade de parafina depositada, deve-se resolver um problema termodinâmico (ESL) para cada um dos nós e tempos utilizados na discretização do domínio. Como resultado, tal problema necessitará ser resolvido milhares de vezes, de acordo com o número de nós na malha em cada simulação realizada. Por essa razão, além de acurado, o modelo empregado para determinar o ESL deverá ser eficiente computacionalmente, isto é, não poderá demandar um tempo computacional muito grande ou o tempo total de execução de cada simulação poderá vir a ser proibitivo. Pelas razões apontadas anteriormente, faz-se necessário realizar uma investigação prévia com o objetivo de determinar as condições necessárias para que se tenha um desempenho otimizado do cálculo das equações do modelo multisólido ¹⁰. Além disso, o elevado número de componentes dos óleos parafínicos requer um estudo mais detalhado sobre aspectos como a robustez (frente a diferentes estimativas iniciais) e tempo de computação . Tais aspectos são importantes no contexto e, até onde os autores puderam determinar, tal análise ainda não foi feita de forma sistemática no problema de deposição de parafinas. Neste trabalho, foi apresentada e discutida a previsão de precipitação de parafinas em óleos parafínicos, considerando o modelo multisólido (múltiplas fases sólidas puras) e duas abordagens distintas para descrever a fugacidade da fase líquida: via equação de estado de Peng-Robinson e via modelos de coeficiente de atividade. Para resolver o SENL, foi utilizado do método de Newton.

O presente trabalho está organizado da seguinte maneira: na seção 2, apresentam-se o modelo multisólido, o teste de estabilidade termodinâmica e as equações de balanço material a serem resolvidas. Na seção 3 descreve-se a modelagem termodinâmica utilizada para este estudo. Na seção 4 são apresentados os algoritmos para solução do SENL. Na seção 5 é discutida a validação dos modelos definidos na modelagem termodinâmica. Já na seção 6 são apresentadas e discutidas as simulações numéricas e feitas considerações sobre os modelos analisados. Por fim, na seção 7 são apresentadas as considerações finais do trabalho.

2 MODELO MULTISÓLIDO

O primeiro passo da aplicação do modelo multisólido corresponde à condução de uma análise de estabilidade ¹⁰. Neste trabalho, foram consideradas as duas abordagens existentes: via modelos de coeficiente de atividade e empregando-se a PR-EOS.

Considerando-se que c é a quantidade de componentes de uma mistura, um componente da mistura poderá existir no estado sólido se a condição expressa pela Eq. (2.1), que representa a abordagem pelo coeficiente de atividade, for satisfeita ⁶:

$z_i^L{y}_i{f}_{pur{o}_i}^L\left(T,P\right)-f_{pur{o}_i}^S\left(T,P\right)\ge 0parai=1,...,c$ (2.1)

Na equação Eq. (2.1), para o i-ésimo componente da mistura, $z_i^L$ representa a composição da alimentação (global), T a temperatura, P a pressão, $f_{pur{o}_i}^S,f_{pur{o}_i}^L$ e γ _i são, respectivamente, a fugacidade do sólido e do líquido puro e o coeficiente de atividade (que será obtido pelos modelos termodinâmicos apresentados na próxima seção).

Ao considerar a segunda abordagem, via PR-EOS, tem-se ¹⁰:

$z_i^L{\varphi }_{i}P-f_{pur{o}_i}^S\left(T,P\right)\ge 0parai=1,...,c$ (2.2)

onde ϕ _i é o coeficiente de fugacidade para um compoente i da mistura calculado via PR-EOS.

Dessa forma, tem-se determinado o número de componentes que existem na fase sólida, representado por c _s .

Na próxima etapa, determina-se o problema de equilíbrio sólido-líquido a ser resolvido (naturalmente, para os componentes que coexistem em ambas as fases). A Eq. (2.3) é referente à abordagem pelo coeficiente de atividade ⁶:

$x_i^L{\gamma }_i{f}_{pur{o}_i}^L\left(T,P\right)-f_{pur{o}_i}^S\left(T,P\right)=0parai=c-c_s+1,...,c$ (2.3)

onde $x_i^L$ representa a fração molar de um componente i em fase líquida (a ser determinanada na solução do problema de coexistência sólido-líquido).

A Eq. (2.4) representa a abordagem via equação PR-EOS ¹⁰:

$x_i^L{\varphi }_{i}P-f_{pur{o}_i}^S\left(T,P\right)=0parai=c-c_s+1,...,c$ (2.4)

A equação de balanço material para os componentes que não se precipitam é a seguinte ¹⁰:

$z_i^L-x_i^L\left[1-\sum _{j=c-c_s+1}^c\frac{n_j^s}{F}\right]=0parai=1,...,c-c_s$ (2.5)

onde $F=1$ é vazão molar da alimentação (base) e $n_j^s$ é o número de mols da fase sólida, enquanto que a equação do balanço material para os componentes que precipitam é ¹⁰:

$z_i^L-x_i^L\left[1-\sum _{j=c-c_s+1}^c\frac{n_j^s}{F}\right]-\frac{n_i^s}{F}=0parai=c-c_s+1,...,c$ (2.6)

A dimensão do SENL em questão é, portanto $c+c_s$. Nazar et al.¹³ recomendam a seguinte expressão para estimar a fugacidade de sólido puro $f_{puroi}^S$:

$$\begin{gathered} \ln \left[\frac{f_{pur{o}_i}^S}{f_{pur{o}_i}^L}\right]=\frac{0,3033}{R}{\mathrm{mml:m}}_i\left(\frac{T_i^f}{T}-1-\ln \left(\frac{T_i^f}{T}\right)\right)-\frac{\Delta H_i^f}{RT}\left(1-\frac{T}{T_i^f}\right)- \\ \frac{\Delta H_i^{tr}}{RT}\left(1-\frac{T}{T_i^{tr}}\right)-\frac{0,0004635}{2R}{\mathrm{mml:m}}_i\left(\frac{{\left(T_i^f\right)}^2}{T}+T-2T_i^f\right) \end{gathered}$$ (2.7)

onde R constante universal dos gases, $\Delta H_i^f$ é a variação da entalpia de fusão, $T_i^f$ é a temperatura de fusão, M _i é a massa molar em (g/mol), $\Delta H_i^{tr}$ é a variação da entalpia de transição e $T_i^{tr}$ é a temperatura transição. A temperatura de transição refere-se a uma transição sólido-sólido que ocorre antes da fusão (e é inferior à temperatura de fusão), sendo dada por ¹⁴:

$T_i^{tr}=366,39775+0,03609{\mathrm{mml:m}}_i-\frac{20879}{{\mathrm{mml:m}}_i}$ (2.8)

A temperatura de fusão para parafinas é dada por ¹⁸:

$T_i^f=374,5+0,02617{\mathrm{mml:m}}_i-\frac{20172}{{\mathrm{mml:m}}_i}$ (2.9)

A variação da entalpia de transição é dada pela expressão proposta por ¹⁴:

$\Delta H_i^{tr}=0,0577{\mathrm{mml:m}}_i{T}_i^f$ (2.10)

Diversas outras expressões para calcular a variação da entalpia de fusão foram abordadas e combinadas e estão detalhadas em Blass ¹.

3 MODELAGEM TERMODINÂMICA

Nesta seção , foram detalhados os modelos para a fase líquida, quais sejam: a equação de estado de Peng-Robinson ¹⁵, o modelo de solução ideal ¹⁷, o modelo UNIFAC ⁹, o modelo UNIQUAC Preditivo ⁶, e o modelo de Wilson Preditivo ¹⁶. De modo a não sobrecarregar a notação , as frações molares na fase líquida nesta seção são apresentadas sem o índice L.

3.1 Equação de Peng-Robinson

Uma das maneiras para descrever a fugacidade de um componente da mistura na fase líquida é empregar uma equação de estado cúbica, como a equação de Peng-Robinson. A expressão explícita na pressão da equação de Peng-Robinson é dada por ¹⁵:

$P=\frac{RT}{V-b}-\frac{a\alpha \left(\omega ,T_r\right)}{V^2+2bV-b^2}$ (3.1)

onde R é a constante universal dos gases, P é a pressão, T é a temperatura absoluta, V é o volume molar da fase. Os parâmetros a e b representam o parâmetro atrativo e o covolume, respectivamente (calculados a partir das propriedades de componentes puros). O parâmetro α é função do fator acêntrico ω e da temperatura reduzida T _r¹⁵.

Da solução da Eq. (3.1) para o volume (com temperatura e pressão especificadas) é possível obter Z, que é o fator de compressibilidade usado na Eq. (3.2), $Z=\frac{PV}{RT}$. Para calcular o coeficiente de fugacidade da fase líquida, emprega-se a expressão ¹⁷:

$$\begin{gathered} \ln {\varphi }_i\left(T,P,x_i^L\right)=\frac{b_i}{b}\left(Z-1\right)-\ln \left(Z-B\right)- \\ \frac{A}{2\sqrt{2}B}\left[\frac{2{\displaystyle \sum _{j=1}^c}{x}_j{a}_{ij}}{a}-\frac{b_i}{b}\right]\ln \left(\frac{Z+2,414B}{Z-0,414B}\right) \end{gathered}$$ (3.2)

em que o coeficiente de fugacidade ϕ _i da Eq. (2.4) é obtido a partir da Eq. (3.2). As definições dos parâmetros a _ij , b _i , A e B podem ser encontradas em Walas ¹⁷.

O coeficiente de fugacidade ϕ _puro,i (T, P) para cada componente puro é dado pela equação PR-EOS ¹⁷:

$\ln {\phi }_{puro,i}\left(T,P\right)=\left(Z-1\right)-\ln \left(Z-B\right)-\frac{A}{2\sqrt{2}B}\ln \left(\frac{Z+2,414B}{Z-0,414B}\right)$ (3.3)

A fugacidade de um componente puro $f_{puroi}^L$ é então calculada por ¹⁷:

$f_{puroi}^L\left(T,P\right)=P{\phi }_{puro,i}\left(T,P\right)$ (3.4)

3.2 Modelo UNIFAC

Para misturas contendo apenas alcanos, pode-se utilizar o modelo UNIFAC ⁹. A expressão para o coeficiente de atividade para uma substância k é:

$\ln {\gamma }_k=1+\ln \left(\frac{{\Phi }_k}{x_k}\right)-\frac{{\Phi }_k}{x_k}-\frac{z}{2}{q}_k\left[1+\ln \left(\frac{{\Phi }_k}{{\theta }_k}\right)-\frac{{\Phi }_k}{{\theta }_k}\right]$ (3.5)

em que o parâmetro z representa o número de coordenação que depende da estrutura molecular, é determinado pela literatura com valor igual a 6 ⁶, θ _k é a fração de área superficial molecular e Φ_k é a fração de volume molecular, os quais são obtidos a partir das seguintes correlações ⁹:

${\theta }_k=x_k{q}_k{\left(\sum _{k=1}^c{x}_k{q}_k\right)}^{-1}{\Phi }_k=x_k{r}_k{\left(\sum _{k=1}^c{x}_k{r}_k\right)}^{-1}$ (3.6)

Para alcanos, r _k é o parâmetro de tamanho (volume) molecular e q _k é o parâmetro molecular de superfície, os quais foram obtidos a partir de Esmaeilzadeh et al.⁶:

$r_k=0,6744C_{nk}+0,4534$ (3.7)

$q_k=0,54C_{nk}+0,616$ (3.8)

onde C _n representa o número de carbonos.

3.3 Modelo de Wilson Preditivo

Neste trabalho, também foi analisado o Modelo de Wilson Preditivo, conforme descrito por Stenby e Coutinho ¹⁶. O coeficiente de atividade para uma espécie k é dado por ⁶:

$\ln {\gamma }_k=1-\ln \left(\sum _{j=1}^c{x}_j{\Lambda }_{kj}\right)-\sum _{i=1}^c\left[\frac{x_i{\Lambda }_{ik}}{{\displaystyle \sum _{j=1}^c}{x}_j{\Lambda }_{ij}}\right]$ (3.9)

Os parâmetros Λ_ki e Λ_ik estão relacionados com as diferenças das energias de interações moleculares e podem ser expressos da seguinte maneira ⁶:

$\ln {\Lambda }_{kj}=-\frac{{\lambda }_{kj}-{\lambda }_{kk}}{RT}$ (3.10)

A energia de interação λ _kk entre duas moléculas idênticas pode ser estimada pela entalpia de sublimação do cristal da molécula pura ⁶:

${\lambda }_{kk}=-\frac{2}{z}\left(\Delta H_k^{sub}-RT\right)$ (3.11)

em que o valor da variação do calor de sublimação $\Delta H_k^{sub}$ foi calculado pela correlação (conhecida como PERT2) dada por Morgan e Kobayashi ¹² e a energia de interação entre as duas moléculas não-idênticas é dada pela seguinte equação ⁶:

${\lambda }_{kj}={\lambda }_{jk}={\lambda }_{jj}$ (3.12)

em que j representa o alcano com a menor cadeia do par kj.

3.4 Modelo UNIQUAC Preditivo

Neste trabalho, foi usado o modelo UNIQUAC Preditivo, desenvolvido por Coutinho ³. O coeficiente de atividade para uma espécie k é dado por ⁶:

$$\begin{gathered} \ln {\gamma }_k=1+\ln \left(\frac{{\Phi }_k}{x_k}\right)-\frac{{\Phi }_k}{x_k}-\frac{z}{2}{q}_k\left[1+\ln \left(\frac{{\Phi }_k}{x_k}\right)-\frac{{\Phi }_k}{x_k}\right]+q_k- \\ {q}_k\ln \left(\sum _{j=1}^c{\theta }_j{\tau }_{jk}\right)-q_k\sum _{j=1}^c\frac{{\theta }_j{\tau }_{kj}}{{\displaystyle \sum _{l=1}^c}{\theta }_l{\tau }_{lj}} \end{gathered}$$ (3.13)

O parâmetro τ _jk é calculado empregando-se a seguinte expressão ⁶:

$\ln {\tau }_{jk}=\frac{{\lambda }_{jk}-{\lambda }_{kk}}{q_{k}RT}$ (3.14)

As quantidades τ _jk e τ _kj foram calculadas de forma similar ao que foi feito para o modelo de Wilson Preditivo. As correlações para os valores de r _k e q _k de n-alcanos são ⁶:

$r_k=0,01148C_{nk}+0,00996$ (3.15)

$q_k=0,0185C_{nk}+0,0211$ (3.16)

3.5 Modelo de Solução Ideal

O modelo de solução ideal assume a fase líquida ideal e, consequentemente, admite-se ${\gamma }_i=1$ para todos os componentes na Eq. (2.3) ¹⁷.

Os resultados indicados como UN/Hi referem-se ao modelo Híbrido, onde o teste de estabilidade termodinâmica conduzido considerando o modelo de solução ideal teve o mesmo número de componentes que tendiam a estar na fase sólida comparado ao modelo UNIFAC. O SENL é gerado usando o modelo UNIFAC; sendo assim, os resultados serão idênticos aos do modelo UNIFAC. Pretende-se, com isso, diminuir o custo computacional ¹.

4 ALGORITMO PARA A SOLUÇÃO DO PROBLEMA NÃO LINEAR

Neste trabalho, o enfoque é exclusivamente investigar a eficiência computacional dos modelos já descritos. Apesar de ser inegavelmente importante, não trataremos aqui das diferentes estratégias numéricas para se resolver a não linearidade presente nos modelos. Tais aspectos estarão presentes em trabalhos futuros. O SENL a ser resolvido tem $c+c_s$ equações e incógnitas. As incógnitas são: $x_i^L$ (c incógnitas) e os $n_i^S$ (c _s incógnitas). Pode-se observar que c _s pode mudar de valor durante o processo de determinação da solução . As equações são dadas pelas Eqs. (2.3) a (2.6), dependendo da abordagem usada, como descrito no Algoritmo 1.

As equações resultantes do equilíbrio sólido-líquido via modelo multisólido, especificadas na seção 2, terão duas abordagens: se for considerada a abordagem pelo coeficiente de atividade, o SENL resultante será dado pela Eq. (2.3) e pelas Eqs. (2.5) e (2.6) representando os balanços materiais. Caso a abordagem escolhida seja via PR-EOS o SENL, as equações serão a Eq. (2.4), referente ao equilíbrio de fases, e as Eqs. (2.5) e (2.6), dos balanços materiais.

4.1 Método de Newton

O procedimento usado (na forma de uma relação de recorrência), até que se atinja uma tolerância definida a priori, é dado pela Eq. (4.1) ⁸:

$\left[J\left(x^k\right)\right]\Delta x=-F\left(x^k\right),k=0,1,2,3,...,n$ (4.1)

em que n é o número máximo de iterações, $\Delta x=x^{k+1}-x^k$, x ^k é o vetor com a estimativa inicial para a composição da fase líquida $\left(x_i^L\right)$ concatenado com o vetor com os componentes existentes na fase sólida $\left(n_i^S\right),J\left(x^{\left(k\right)}\right)$ é a matriz Jacobiana obtida das equações do modelo multisólido, $F\left(x^k\right)$ é o conjunto de equações do SENL resultante e k é o contador de iterações. Os algoritmos para a solução do problema são apresentados a seguir.

Algoritmo 1
Algoritmo para escolha da abordagem e teste de estabilidade

Dados de Entrada:T , P, c, $z_i^L$ com $i=1,...,c$ e uma estimativa inicial de $x_i^L$
1:	Calcular $f_{puroi}^L,parai=1,...,c$, via Eq. (3.4)
2:	Calcular $f_{puroi}^S,parai=1,...,c$, via Eq. (2.7)
3:	for all $i=1,...,c$ do
4:		if (Abordagem via coeficiente de atividade para calcular $f_i^L\left(T,P,z_l^L\right)$) then
5:			Escolher: UNIFAC, UNIQUAC preditivo, Wilson preditivo ou Solução Ideal
6:		else
7:			PR-EOS
8:		end if
9:	end for
10:	for all $i=1,...,c$ do
11:		if $f_i^L\left(T,P,z_i^L\right)-f_{pur{o}_i}^S\ge 0$ then
12:			A i-ésima componente da mistura existe no estado sólido
13:		else
14:			A mistura existe no estado líquido
15:		end if
16:	end for

O procedimento iterativo é interrompido quando um critério de parada estabelecido a priori é satisfeito (o erro absoluto inferior a uma dada tolerância no presente trabalho). A tolerância para interrupção do algoritmo numérico usada em todos os casos é $\epsilon ={10}^{-5}$. A pressão é mantida constante, não tendo influência significativa nos resultados. O Algoritmo 1 tem como objetivo escolher uma das abordagens para calcular a fugacidade do componente i na fase líquida, bem como efetuar o teste de estabilidade termodinâmica. O Algoritmo 2 apresenta a construção do SENL, de acordo com a abordagem (modelo), juntamente com os balanços materiais e o método numérico para resolução do sistema de equações.

Como apontado na introdução deste trabalho, há a necessidade de se ter um modelo acurado para determinar o ESL. Além de preciso, tal modelo tem que ser eficiente do ponto de vista computacional. Dois aspectos têm papel relevante na determinação de tal eficiência. O primeiro é o modelo propriamente dito, pois modelos diferentes levam tempos diferentes para alcançar soluções semelhantes, como ficará demonstrado na seção de resultados. O segundo aspecto é como resolver o sistema de equações não lineares existentes nesses modelos. Por exemplo, se o caminho escolhido for usar o método de Newton, dependendo da estratégia numérica empregada na solução do sistema de equações (se forem ou não empregados precondicionadores, sobre qual entre os inúmeros algoritmos iterativos existentes são usados para resolver o sistema de equações lineares etc.), o tempo gasto para determinar o ESL também irá variar. Não há dúvidas de que a variação será significativa, como já observado por Guterres ⁷, que identificou significantes variações do tempo de CPU, dependendo da estratégia computacional escolhida para resolver o sistema de equações não-lineares da equação de Richards. Resumindo, o método para determinar o ESL mais adequado deverá ter um algoritmo eficiente, combinado com técnicas numéricas sofisticadas para resolver o sistema de equações não lineares.

Algoritmo 2
Algoritmo que define o SENL e o Método de solução

	for all $i=1,...,c$ do
2:		if (Usar da abordagem: Eq. (2.3) -- Eq. (2.4) then
			SENL para $i=c-c_s+1,...,c$
4:			Calcular o balanço para os componentes que não se precipitam Eq. (2.5)
			SENL para $i=1,...,c-c_s$
6:			Calcular o balanço para os componentes que precipitam Eq. (2.6)
			SENL para $i=1,...,c-c_s$
8:		end if
	end for
10:	Para resolver o SENL de $c+c_s$ equações, usar o Método de Newton, via Eq. (4.1)
	while $\left(\mathrm{mml:m}ax\left|x_i^{k+1}-x_i^k\right|>0,00001i=1,...,c+c_s\right)$ do
12:		Aplicar a relação de recorrência do Método de Newton, Eq. (4.1)
	end while

Pelas razões supracitadas, foi escolhido um método direto, a eliminação de Gauss, para resolver o sistema de equações lineares associado ao Método de Newton que está representado pela Eq. (4.1). Nota-se que não há dúvidas de que ao se resolver o sistema linear dessa forma, opta-se por um método que não é o mais eficiente. Contudo, exatamente por ser um método direto e não um método iterativo, o tempo de CPU gasto para resolver um sistema de equações lineares dependerá basicamente do número de equações do sistema. Dessa forma, ao se comparar o tempo para se alcançar a solução de uma determinada mistura modelada com duas formulações distintas, a diferença de tempo de CPU para se alcançar a solução do problema ESL deverá ser atribuída somente aos modelos propriamente ditos e não em função de que, por alguma razão, um determinado modelo resolveu as equações do sistema iterativamente em menor tempo computacional.

Os elementos da matriz jacobiana presentes na Eq. (4.1) foram calculados utilizando uma simples expressão de diferenças finitas:

$J_{ij}=\frac{F_i\left(x_1,x_2,...,x_j+\Delta x_j,...,x_c\right)-F_i\left(x_1,x_2,...,x_c\right)}{\Delta x_j}$ (4.2)

5 VALIDAÇÃO DOS CÓDIGOS

De modo a apresentar uma validação do código computacional, calculou-se o percentual de massa na fase sólida para diversas misturas onde existem dados experimentais disponíveis. A pressão considerada foi de 1 atm (foi verificado que não há influência considerável para esse tipo de problema, como esperado, mesmo quando se empregou a equação de estado). Os resultados numéricos foram comparados com dados experimentais extraídos da literatura ⁵. Na Tabela 1 estão apresentadas as composições mássicas das misturas que foram utilizadas para a validação dos códigos. As misturas têm em sua composição ${\mathrm{mml:m}}_1=17$ componentes, ${\mathrm{mml:m}}_2=15$ componentes, ${\mathrm{mml:m}}_3=11$ componentes e ${\mathrm{mml:m}}_4=7$ componentes ⁵.

Tabela 1:
Composição percentual (base mássica) das misturas analisadas. ⁵

Componentes	M 1	M 2	M 3	M 4
n − C10	64,250	64,650	65,590	66,650
n − C18	4,872	5,383	6,958	10,570
n − C19	4,411	4,874	6,301	9,571
n − C20	3,981	4,400	5,688	8,641
n − C21	3,587	3,963	5,124	−
n − C22	3,221	3,561	4,603	−
n − C23	2,889	3,193	−	−
n − C24	2,586	2,858	−	−
n − C25	2,310	−	−	−
n − C29	1,451	−	−	−
n − C30	1,290	1,424	−	−
n − C31	1,142	1,262	−	−
n − C32	1,011	1,119	1,444	−
n − C33	0,895	0,989	1,280	−
n − C34	0,790	0,873	1,129	1,716
n − C35	0,698	0,771	0,998	1,516
n − C36	0,616	0,682	0,881	1,337

Os dados obtidos pelos modelos estudados foram comparados com os dados experimentais ⁵. Para todas as misturas analisadas neste trabalho foi calculado o desvio absoluto $\left|x_{calculado}-x_{experi\mathrm{mml:m}ental}\right|$, que resulta na diferença entre o percentual de massa sólida predito pelos modelos versus o resultado experimental para cada temperatura. Para a mistura M ₁, a temperatura onde ocorreu o desvio máximo foi de 263,15 K, enquanto que a temperatura onde houve desvio mínimo foi 296,15 K.

A Tabela 2 apresenta os desvios das quantidades de fase sólida obtidos em cada modelo nas temperaturas de 263,15 K e de 296,15 K, para a mistura M ₁.

Tabela 2:
Desvios para a mistura M ₁.

Modelos	Máximo	Mínimo
Solução Ideal	4,13	1,35
EOS-PR	3,70	1,96
UN/Hí	3,37	0,73
UNIQUAC	1,76	0,48
Wilson	3,95	0,36

Como observado na Tabela 2, os modelos que tiveram resultados mais aderentes aos experimentais foram o UNIQUAC (para a temperatura onde houve desvio máximo) e Wilson (para a temperatura de desvio mínimo). Na Figura 1, estão representados os percentuais de parafina na fase sólida, para a mistura M ₁ em função de diferentes temperaturas.

Figura 1:
Porcentual de massa sólida para a M ₁.

A Tabela 3 representa os desvios máximos e mínimos (percentuais de fase sólida) comparados com os experimentais para a mistura M ₂. Observa-se que os resultados mais próximos dos dados experimentais foram: UNIQUAC para a temperatura de desvio máximo e o modelo de Wilson para a temperatura onde houve desvio mínimo.

Tabela 3:
Desvios para a mistura M ₂.

Modelos	Máximo	Mínimo
Solução Ideal	4,46	2,43
EOS-PR	3,54	1,15
UN/Hí	3,81	0,55
UNIQUAC	2,30	1,53
Wilson	4,34	0,22

Na Figura 2 podem-se observar os resultados da percentagem de parafina na fase sólida para a mistura M ₂ utilizando-se dos diversos modelos descritos neste texto. Para a mistura M ₂, a temperatura para a desvio máximo é de 263,05 K e para o desvio mínimo é de 273,35 K.

Figura 2:
Percentual de massa sólida para a M ₂.

Para a mistura M ₃, a temperatura para o desvio máximo foi de 263,05 K e para o desvio mínimo foi de 303,15 K. A Tabela 4 representa os desvios máximos e mínimos, sobre os quais é possível observar que o modelo UNIQUAC foi o que apresentou a menor diferença, para ambas temperaturas.

Tabela 4:
Desvios para a mistura M ₃.

Modelos	Máximo	Mínimo
Solução Ideal	5,29	1,57
EOS-PR	4,24	1,63
UN/Hí	4,60	0,72
UNIQUAC	3,92	0,15
Wilson	5,26	0,40

A Figura 3 representa a porcentagem de parafina na fase sólida para a mistura M ₃, em função da temperatura. Novamente, há concordância entre os resultados experimentais e os simulados.

Figura 3:
Percentual de massa sólida para a M ₃.

Para a mistura M ₄, a temperatura para o desvio máximo é de 273,15 K e o desvio mínimo ocorre em 299,15 K. A Tabela 5 representa os desvios máximos e mínimos comparados com experimentais, para ambas as temperaturas. Nesse caso, o modelo UNIQUAC apresentou melhor resultado para a temperatura onde ocorre a diferença máxima, enquanto o modelo de Wilson apresentou o menor desvio para a temperatura onde há a diferença mínima.

Tabela 5:
Desvios para a mistura M ₄.

Modelos	Máximo	Mínimo
Solução Ideal	4,46	0,45
EOS-PR	2,99	0,47
UN/Hí	2,13	0,34
UNIQUAC	0,68	0,34
Wilson	6,41	0,27

A Figura 4 representa a porcentagem de parafina na fase sólida, como função da temperatura para os componentes da M ₄, onde percebe-se claramente a influência da temperatura na precipitação , pela pequena porcentagem de parafina precipitada em temperaturas altas.

Figura 4:
Percentual de massa sólida para a M ₄.

6 RESULTADOS E DISCUSSÕES

Nesta seção estão detalhados os resultados das simulações dos percentuais de massa sólida em função da temperatura para os diferentes modelos para a fase líquida apresentados na seção 3, usando da combinação Solução Ideal/UNIFAC (denominado modelo Híbrido) para ganho computacional. O modelo identificado no trabalho como híbrido (denominado UN/Hi) foi constituído por duas etapas: (i) a etapa de teste de estabilidade termodinâmica, feita com o modelo de Solução Ideal e (ii) o cálculo de equilíbrio de fases feito via modelo UNIFAC. Com o modelo definido, foram analisados alguns cenários.

No cenário 1, foi avaliada a influência da temperatura no percentual de massa da fase sólida e o tempo de execução para cada modelo estudado. A Tabela 6 apresenta o tempo de execução total de cada modelo, considerando uma variação de temperatura de 5 em 5 K, com uma temperatura inicial de 265 K.

Tabela 6:
Tempo de execução para cada modelo em segundos.

Modelos	Método de Newton
Solução Ideal	1,2
EOS-PR	17,1
UNIFAC	18,1
UNIQUAC	2,9
Híbrido	2,8
Wilson	19,1

A Figura 5 ilustra os percentuais de fase sólida para as diversas temperaturas e modelos. Naturalmente, a composição global é mantida constante (conforme apresentada na Tabela 1). Ainda assim, e de acordo com o que foi verificado na validação dos códigos, os resultados para os diferentes modelos estudados são consistentes. Os resultados indicados como UN/Hi referem-se ao modelo Híbrido. Como tal simplificação foi feita somente na etapa de cálculo de estabilidade termodinâmica (considerando o modelo de Solução Ideal), os resultados serão idênticos aos do modelo UNIFAC. Tal consideração é válida para todas as figuras envolvendo o cálculo de percentuais de massa sólida em função da temperatura. Os resultados são pertinentes para a mistura M ₁, tendo em conta que é a mistura M ₁ a que apresenta maior número de componentes.

Figura 5:
Percentual de massa sólida para a mistura M ₁.

No segundo cenário, foi empregada uma função de comportamento senoidal (para testar o comportamento do arcabouço computacional) como dado de entrada para temperatura e a composição inicial utilizada foi a da mistura M ₁, apresentada na Tabela 1. O intervalo [0; 42], totalizando em quarenta e três pontos, foi definido como um certo conjunto de pontos da “malha”. Como observado na Figura 6, a proposta foi ter um maior número de pontos em relação ao cenário 1 e, como resultado, tem-se uma medida mais clara dos tempos de computação envolvidos em possíveis simulações de escoamentos multifásicos. Os resultados geram comportamentos esperados, ou seja, pode-se verificar na Figura 6 que o comportamento da precipitação de parafina sólida produzido é também oscilatório, resultado da entrada de temperatura senoidal.

Figura 6:
Percentual de massa sólida para a mistura M ₁.

Na Tabela 7, podem-se observar os tempos totais de execução . Considerando a mistura M ₁, observou-se que a utilização do método híbrido produz algum ganho computacional em relação ao modelo UNIFAC.

Tabela 7:
Tempo de execução para cada modelo em segundos.

Modelos	Método de Newton
Solução Ideal	7,2
EOS-PR	87,4
UNIFAC	92,3
UNIQUAC	19,1
Híbrido	16,2
Wilson	48,2

Para o cenário 3, as frações molares das mistura M ₁ da Tabela 1 foram geradas aleatoriamente, em uma distribuição uniforme, resultando em 43 frações molares diferentes para cada componente da Mistura M ₁. Ou seja, no cenário 3, cada ponto tem uma composição diferente e foram obtidos os perfis de quantidades de parafina sólida, variando-se simultaneamente a temperatura e a composição da mistura. Assim, esse cenário é um teste mais severo para a análise de robustez dos códigos e também para uma avaliação mais genérica dos tempos de computação . Os resultados dessa simulação podem ser vistos na Figura 7, que apresenta os dados do percentual de massa na fase sólida em função da temperatura e composição em cada ponto, enquanto a Tabela 8 traz o tempo de execução total de cada método para a mistura M ₁. Verifica-se que o modelo PR-EOS divergiu em pelo menos um ponto. Os resultados justificam uma investigação para o uso de outros métodos para resolver o sistema não linear resultante. Portanto, no contexto do desenvolvimento de simuladores (por exemplo, em problemas de escoamento), o uso de técnicas de globalização de métodos do tipo Newton deve ser testado, de modo a promover um aumento de robustez da ferramenta computacional.

Figura 7:
Percentual de massa sólida para a mistura M ₁.

Tabela 8:
Tempo de execução para cada modelo em segundos.

Modelos	Método de Newton
Solução Ideal	7,0
PR-EOS	Divergiu
UNIFAC	97,1
UNIQUAC	18,8
Híbrido	15,7
Wilson	85,1

7 CONCLUSÕES

Neste trabalho foi apresentada e discutida a previsão de precipitação de ceras em óleos parafínicos, considerando o modelo multisólido (múltiplas fases sólidas puras), com diferentes abordagens para a descrição das fugacidades da fase líquida (PR-EOS, Solução Ideal, modelo de Wilson preditivo, modelo UNIFAC e modelo UNIQUAC preditivo). Uma abordagem híbrida, formada pelo modelo de Solução Ideal para a etapa de teste de estabilidade termodinâmica e pelo modelo UNIFAC na solução do problema de equilíbrio sólido-líquido foi proposta e implementada com sucesso. Assim como a combinação dos modelos, foram resolvidos problemas de equilíbrio para diferentes temperaturas e composições. Os tempos de computação e a robustez do método numérico utilizados para resolver o sistema não linear foram avaliados. Os resultados apresentados e discutidos na seção 6 indicam que: (a) Conforme esperado, a complexidade do modelo termodinâmico que representa a mistura influencia diretamente no tempo de computação e na robustez do método numérico (cenários abordados neste trabalho que não foram explorados por outros autores); (b) Será necessário o uso de outros tipos de métodos numéricos que fazem o uso de técnicas de globalização para a resolução dos sistemas não-lineares resultantes. Como pode ser observado na Tabela 8, para o modelo PR-EOS houve divergência.

Finalmente, um ponto central a ser avaliado quando da implementação de problemas de equilíbrio sólido-líquido é que essas combinações e análises não se encontram na literatura especializada. As análises realizadas apontam para uma promissora possibilidade de construção de arcabouços computacionais robustos e computacionalmente eficientes.

REFERÊNCIAS

• ¹
  L. Blass. “Experimentos Numéricos para Deposição de Parafinas com Modelo Multisólido”. Doutorado em Modelagem Computacional, Universidade do Estado do Rio de Janeiro (2015).
• ²
  J.A.P. Coutinho. Predictive UNIQUAC: A New Model for the Description of Multiphase Solid Liquid Equilibria in Complex Hydrocarbon. Industrial & Engineering Chemistry Research, (3) (1998), 4870- 4875.
• ³
  J.A.P. Coutinho. Predictive local composition models: NRTL and UNIQUAC and their application to model solid-liquid equilibrium of n-alkanes. Fluid Phase Equilibria, 158 (1999), 447-457. doi: 10.1016/S0378-3812(99)00085-0.
  » https://doi.org/10.1016/S0378-3812(99)00085-0
• ⁴
  J.A.P. Coutinho, K. Knudsen, S.I. Andersen & E.H. Stenby. A local composition model for paraffinic solid solutions. Chemical Engineering Science, 51(12) (1996), 3273-3282. doi: 10.1016/0009-2509(95)00397-5. URL http://linkinghub.elsevier.com/retrieve/pii/0009250995003975
  » https://doi.org/10.1016/0009-2509(95)00397-5» http://linkinghub.elsevier.com/retrieve/pii/0009250995003975
• ⁵
  C. Dauphin, J.L. Daridon, J.P. Coutinho, P. Baylère & M. Potin-Gautier. Wax content measurements in partially frozen paraffinic systems. Fluid Phase Equilibria, 161 (1999), 135-151.
• ⁶
  F. Esmaeilzadeh, J.F. Kaljahi & E. Ghanaei. Investigation of different activity coefficient models in thermodynamic modeling of wax precipitation. Fluid Phase Equilibria, 248 (2006), 7-18.
• ⁷
  M.X. Guterres. “Avaliação dos Algoritmos de Picard-Krylov e Newton-Krylov na Solução da Equação de Richards”. Doutorado em Modelagem Computacional, Universidade do Estado do Rio de Janeiro (2013).
• ⁸
  C.T. Kelley. “Solving nonlinear equations with Newton’s method”. SIAM Society of Industrial and Applied Mathematics, Philadelphia (2003).
• ⁹
  B.L. Larsen, P. Rasmussen & A. Fredenslund. A modified UNIFAC group-contribution model for prediction of phase equilibria and heats of mixing. Industrial & Engineering Chemistry Research, 26(3) (1987), 2274-2286.
• ¹⁰
  C. Lira-Galeana, A. Firoozabadi & J.M. Prausnitz. Thermodynamics of wax precipitation in petroleum mixtures. AIChE Journal, 42(1) (1996), 239-248. doi: 10.1002/aic.690420120. URL http://doi.wiley.com/10.1002/aic.690420120
  » https://doi.org/10.1002/aic.690420120» http://doi.wiley.com/10.1002/aic.690420120
• ¹¹
  R. Monroy-Loperena. An efficient method to calculate multisolid-(wax) precipitation. Fluid Phase Equilibria, 348 (2013), 70-78. doi: 10.1016/j.fluid.2013.03.005. URL http://linkinghub.elsevier.com/retrieve/pii/S0378381213001222
  » https://doi.org/10.1016/j.fluid.2013.03.005» http://linkinghub.elsevier.com/retrieve/pii/S0378381213001222
• ¹²
  D.L. Morgan & R. Kobayashi. Extension of Pitzer CSP models for vapor pressures and heats of vaporization to long-chain hydrocarbons. Fluid Phase Equilibria, 94 (1994), 51-87.
• ¹³
  A.R.S. Nazar, B. Dabir & M.R. Islam. A Multi-Solid Phase Thermodynamic Model for Predicting Wax Precipitation in Petroleum Mixtures. Energy Sources, 27(1-2) (2005), 173-184. doi: 10.1080/00908310490448253. URL http://dx.doi.org/10.1080/00908310490448253
  » https://doi.org/10.1080/00908310490448253» http://dx.doi.org/10.1080/00908310490448253
• ¹⁴
  D. Nichita, L. Goual & A. Firoozabadi. Wax Precipitation in Gas Condensate Mixtures (2001). doi: 10. 2118/74686-PA.
  » https://doi.org/10. 2118/74686-PA
• ¹⁵
  D.Y. Peng & D.B. Robinson. A new two-constant equation of state. Industrial & Engineering Chemistry Fundamentals, 15(1) (1976), 59-64. doi: 10.1021/i160057a011.
  » https://doi.org/10.1021/i160057a011
• ¹⁶
  E.H. Stenby & J.A.P. Coutinho. Predictive Local Composition Models for Solid/Liquid Equilibrium in n-Alkane Systems: Wilson Equation for Multicomponent Systems. Industrial & Engineering Chemistry Research, 35(3) (1996), 918-925.
• ¹⁷
  S.M. Walas. “Phase Equilibria in Chemical Engineering”. Butterworth-Heinemann, 1st ed. (1985), 688 pp.
• ¹⁸
  K. Won. Thermodynamic calculation of cloud point temperatures and wax phase compositions of refined hydrocarbon mixtures. Fluid Phase Equilibria, 53 (1989), 377-396. doi: 10.1016/0378-3812(89)80104-9.
  » https://doi.org/10.1016/0378-3812(89)80104-9
• ¹⁹
  K.W. Won. Thermodynamics for solid solution-liquid-vapor equilibria: wax phase formation from heavy hydrocarbon mixtures. Fluid Phase Equilibria, 30 (1986), 265-279. doi: 10.1016/0378-3812(86)80061-9. URL http://dx.doi.org/10.1016/0378-3812(86)80061-9
  » https://doi.org/10.1016/0378-3812(86)80061-9» http://dx.doi.org/10.1016/0378-3812(86)80061-9

Datas de Publicação

• Publicação nesta coleção
  30 Abr 2020
• Data do Fascículo
  Jan-Apr 2020

Histórico

• Recebido
  25 Ago 2018
• Aceito
  16 Jun 2019