Simulação Numérica de Escoamento Eletroosmótico Usando o Modelo Constitutivo de Phan-Thien-Tanner

RESUMO

Neste trabalho será investigado o comportamento de escoamentos de fluidos newtonianos e não-newtonianos em microcanais. O problema não-newtoniano, consiste em resolver as equações que regem o movimento para o caso de um escoamento de fluidos cujas propriedades reológicas possam ser estudadas pelo modelo constitutivo de Phan-Thien-Tanner, como por exemplo os materiais poliméricos. Uma das características interessantes de alguns destes materiais é que eles podem ser misturados com solventes apropriados, como uma solução eletrolítica, e o resultado é que este fluido como um todo passa a ter propriedades elétricas. Assim, além das propriedades viscoelásticas, será investigada a eletrocinética do escoamento, que é diretamente influenciado pela aplicação de um campo elétrico externo. Em particular o fenômeno de eletroosmose será estudado por meio de simulações numéricas em canais planos. O movimento das cargas na solução é descrito pelas equações de Poisson-Nernst-Planck e para resolver numericamente este problema será aplicado o método das diferenças finitas generalizadas. O código para as simulações de escoamentos eletroosmóticas foi implementado como uma parte do sistema chamado HiG-Fow.

Palavras-chave:
HiG-Flow; viscoelástico; eletroosmótico

ABSTRACT

In this work the behavior of newtonian and non-Newtonian fluids in microchannels will be investigated. The non-Newtonian problem consists in solving the governing equations of the movement take into account a fluid flow whose rheological properties can be studied by the Phan-Thien-Tanner constitutive model, for example the polymeric materials. One of the interesting features of these materials is that they can be mixed with appropriate solvents, such as an electrolyte solution, and the resulting fluid has electrical properties. Thus, in addition to the viscoelastic properties, the electrokinetics properties of the flow will be investigated, which is directly influenced by the application of an external electric field. In particular the phenomenon of electro-osmosis will be studied through numerical simulations in flat channels. The motion of the charges in the solution is described by the Poisson-Nernst-Planck equations and to solve numerically this problem will be applied the generalized finite difference method. The code for the simulations of electroosmotic flows was implemented as a part of the HiG-Fow system.

Keywords:
HiG-Flow; viscoelastic; electro-osmotic; finite differences

1 INTRODUÇÃO

Há mais de 200 anos, é conhecido que determinadas soluções líquidas podem escoar através de um pequeno canal pela aplicação de uma diferença de potencial elétrico em suas extremidades. Esta área de estudo, ou técnica de manipulação do fluido é chamada de Eletrocinética. Neste âmbito, existem várias aplicações em áreas de pesquisa como química, engenharia e biomédica ^{13), (9}. A superfície em contato com uma solução eletrolítica desempenha um papel importante na eletrocinética de dispositivos microfluidicos. Vórtices formam-se em torno de uma superfície condutora que está em contato com a solução aquosa quando um campo elétrico externo é aplicado. Estes vórtices causam perturbações do fluxo no microcanal, alterando assim o campo de velocidades do fluido. A eletrohidrodinâmica é a conexão entre as teorias de eletromagnetismo e hidrodinâmica. Assim, a dinâmica de fluidos, eletrocinética e eletroquímica formam a base para o estudo da eletrohidrodinâmica. Na teoria eletromagnética, as equações de Maxwell são bem conhecidas. Uma revisão geral e maior detalhes sobre princícios físicos dos efeitos eletrocinéticos pode ser encontrada em ⁶.

O fenômeno de eletroosmose foi primeiramente demonstrado por Reuss em 1809 ²⁷. Reuss demonstrou que o fluxo de fluido polar é conduzido pelo campo elétrico externo aplicado entre a entrada e saída de um canal. Ele também notou que o campo elétrico agia diretamente sobre os íons existentes próximo as paredes do canal. Posteriormente Burgreen e Nakache ⁷ realizaram estudos em capilares ultrafinos, observando o efeito do potencial de superfície sobre o escoamento. Um ano depois, em 1965 Rice e Whitehead publicaram o artigo ²⁸ sobre os mesmos estudos, agora em capilares cilíndricos. Em 1879, Helmholtz ²⁰ relatou os parâmetros elétricos e do escoamento para o transporte eletrocinético e deu o nome de teoria da dupla camada elétrica (Eletrical Double layer-EDL). Posteriormente, Goy ¹⁹ e Chapman ¹⁰ contribuiram nos estudos da distribuição de cargas no fluido próximo as paredes do capilar. Em 1918 von Smoluchowski ²⁹ realizou estudos em capilares, generalizando a teoria proposta por Helmholtz. Outra grande contribuição foi dada por Debye e Hückel em 1923 ¹¹. Eles determinaram a concentração iônica na solução por meio da linearização da distribuição de Boltzmann para a energia. Posteriormente varios estudos relacionados a efeitos eletrocinéticos em fluidos newtonianos foram realizados. Burgreen e Nakache ⁷ realizaram estudos em capilares ultrafinos, observando o efeito do potencial de superfície sobre o escoamento. Um ano depois, em 1965 Rice e Whitehead publicaram o artigo ²⁸ sobre os mesmos estudos, agora em capilares cilíndricos. Análises e simulações numéricas de escoamentos eletroosmóticos começaram a ser publicadas por Yang e Li ³⁵ e Patankar e Hu ²³. O método de diferenças finitas foi aplicado por Ermakov et al. para resolver escoamentos eletroosmóticos em duas dimensões para fluidos newtonianos. No ano de 2000, Bianchi et al ⁵ resolveram o problema de escoamento eletroosmótico em uma junção-T usando o método de elementos finitos. Soluções analíticas para escoamentos newtonianos em canais foram demonstradas por Dutta e Beskok ¹⁴. Em 2002 Lin et al. ²¹ resolveram numericamente a equação de Nernst-Planck para o transporte iônico junto as equações de Navier-Stokes para fluidos newtonianos. O estudo de escoamento eletroosmótico para fluidos não-newtonianos é mais recente, devido a dificuldade imposta dependendo do modelo constituvo a ser utilizado. Em 2008 Park e Lee ²² determiram a velocidade eletroosmótica para escoamentos viscoelásticos enquanto Zhao et al. ³⁶ realizavam estudos de escoamentos eletroosmóticos de fluidos modelados pela lei de potência. Tang et al. ³³ também realizaram estudo numérico de fluidos não-newtonianos em escoamentos eletroosmóticos. Afonso et al. ¹ demonstraram a solução analítica para escoamentos viscoelásticos em um canal de seção circular, levando em conta um gradiente de pressão não nulo e o efeito eletrocinético, onde um dos modelos constitutivo usado foi o de Phan-Thien/Thanner simplificado, sPTT. Em 2015, Peng et al. ²⁴ estudaram os efeitos de um gradiente de concentração dentro de um canal, levando em conta uma mistura de diferentes soluções eletrolíticas. Recentemente, Song et al. ³⁰ realizaram estudos de instabilidades numéricas em mistura de uma solução ferrosa com agua em um canal plano.

Um dos problemas numéricos que surge quando estudamos escoamentos eletroosmóticos consiste essencialmente em descrever a concentração de cargas em função do tempo, isto é, precisamos obter a densidade de cargas elétricas no fluido. Esta densidade de cargas é dada pela solução da equação de Nernst-Planck, que é uma equação de transporte dos íons. Além disso, esta equação é acoplada ao potencial induzido nas paredes do canal pelo qual ocorre o escoamento. Assim, temos que resolver uma equação de advecção-difusão, problema similar ao que aparece nas equações de Navier-Stokes. Para resolver este problema numericamente, o domínio computacional do escoamento é dado por uma malha hierarquica gerada através da Hig-tree ³¹, por exemplo, no caso bidimensional é uma quad-tree generalizada ¹⁷. Em geral, as malhas hierárquicas impõem dificuldades no esquema numérico com base em aproximações cartesianas e requer o uso de interpolações espaciais em pontos desconhecidos do estêncil. As interpolações das propriedades no centro das faces e no centro das células são feitas pela técnica de mínimos quadrados móveis, que usa um dado conjunto de pontos onde a propriedade é conhecida para estimar um valor desconhecido em um ponto vizinho. As equações de Navier-Stokes juntamente com o modelo constitutivo de Phan-Thien-Tanner (PTT) são discretizadas pelo método das diferenças finitas generalizadas usando o sistema HiG-Flow ⁸, que permite que ”compartimentos de código”viscoelástico, newtoniano, eletroosmótico, sejam criados e depois usados juntos para a simulação numérica. Os sistemas lineares resultantes das discretização são resolvidos por solvers usando a biblioteca PETSc (Portable, Extensible Toolkit for Scientific Computation) ⁴.

2 DESCRIÇÃO DO PROBLEMA

Como foram estudados tanto escoamentos de fluidos newtonianos quanto viscoelásticos, iremos apresentar as equações governantes de forma geral para o caso viscoelástico e, veremos que o caso newtoniano é descrito por uma simplificação dessas equações. Além disso, as equações serão apresentadas na forma adimensional, pois, o código está implementado desta forma. Neste sentido, os parâmetros, operadores e propriedades adimensionalizados serão identificados por um superescrito asterisco ”∗”. Consideramos que o escoamento é incompressível, laminar, isotérmico e com viscosidade constante. Em particular, escoamentos em microfluidica tem baixas velocidades, isto é, bem menores que a velocidade do som, por isso é comum considerá-lo como sendo incompressível e laminar, sem comprometer a acurácia do método de solução. Assim, as equações que regem o escoamento na forma dimensional são:

$\nabla ·\mathbf{u}=0,$ (2.1)

$\rho \frac{\partial \mathbf{u}}{\partial t}+\rho \mathbf{u}·\nabla \mathbf{u}=-\nabla p+{\eta }_0{\nabla }^2\mathbf{u}+\nabla ·\mathbf{S}+\mathbf{F},$ (2.2)

$\mathbf{T}=2{\eta }_p\mathbf{D}+\mathbf{S},$ (2.3)

O sistema HiG-Flow fora implementado de maneira que todas as equações são resolvidas na forma adimensional. Para adimensionalizar as equações governantes (2.1), (2.2) e (2.3) considera-se as grandezas adimensionais

$p^{*}=\frac{p}{\rho U^2}{u}^{*}=\frac{u}{U}{v}^{*}=\frac{v}{U}{x}^{*}=\frac{x}{H}{t}^{*}=\frac{tU}{H}{S}^{*}=\frac{S}{\rho U^2}{T}^{*}=\frac{T}{\rho U^2},$

onde ρ, H e U são valores dimensionais de referência para densidade, comprimento e velocidade respectivamente. Neste sentido, ${\nabla }^{*}=H\nabla$ e logo ${\mathbf{D}}^{*}=H/U\mathbf{D}$ permitindo escrever as equações na forma adimensional:

${\nabla }^{*}·{\mathbf{u}}^{*}=0,$ (2.4)

$\frac{\mathit{\partial }{\mathbf{u}}^{*}}{\mathit{\partial }{t}^{*}}+{\mathbf{u}}^{*}·{\nabla }^{*}{\mathbf{u}}^{*}=-{\nabla }^{*}{p}^{*}+\frac{1}{Re}{\nabla }^{2*}{\mathbf{u}}^{*}+{\nabla }^{*}·{\mathbf{S}}^{*}+{\mathbf{F}}^{*},$ (2.5)

${\mathbf{T}}^{*}=\frac{2(1-\beta )}{Re}{\mathbf{D}}^{*}+{\mathbf{S}}^{*},$ (2.6)

onde u é o campo vetorial de velocidade, t é o tempo, p é a pressão hidrodinâmica, $Re=\rho UH/{\eta }_0$ é o número de Reynolds, U é a velocidade característica do fluxo, H é a largura característica do canal por onde ocorre o escoamento, ρ a densidade de massa e η ₀ denota a viscosidade dinâmica total ${\eta }_0={\eta }_s+{\eta }_p$. A proporção de solvente é dada por $\beta =\frac{{\eta }_s}{{\eta }_0}$. O tensor taxa de deformação é denotado por $\mathbf{D}=\frac{1}{2}\left(\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^T\right)$, T o tensor das tensões elástica, também chamado neste texto de tensor polimérico e F é o termo relativo as forças externas aplicadas ao fluido, o qual pode depender do campo gravitacional, de um campo elétrico ou outros campos de forças. A contribuição polimérica T é definida em termos de um novo tensor, S, na equação 2.3. Esta relação entre T e S permite que o termo difusivo ${\eta }_o{\nabla }^2\mathbf{u}$ apareça na equação de movimento. Caso o escoamento seja do tipo viscoelástico, este termo difusivo some, sendo computada a contribuição do solvente se houver, e a contribuição polimérica. Já se o escoamento for do tipo newtoniano, o tensor S é nulo, consequentemente ${\eta }_o={\eta }_s$ e as equações 2.2, 2.3 são consistentes. A equação de movimento (2.5) é um balanço de forças que agem sobre o fluido. A evolução no tempo do tensor polimérico é dada por:

$\frac{{\partial }^{*}{\mathbf{T}}^{*}}{\partial t^{*}}+\left({\mathbf{u}}^{*}·{\nabla }^{*}\right){\mathbf{T}}^{*}-\left[{\left({\nabla }^{*}{\mathbf{u}}^{*}\right)}^T·{\mathbf{T}}^{*}+{\mathbf{T}}^{*}·{\nabla }^{*}{\mathbf{u}}^{*}\right]=\frac{1}{De}{\mathbf{M}}^{*}\left({\mathbf{T}}^{*}\right),$ (2.7)

onde $De=\lambda U/H$ é o número de Deborah e λ é o tempo de relaxação do fluido. No contexto geral do código, a função M * (T *) é genérica, podendo ser definida mais tarde. O sistema HiG-Flow carrega a função apropriada, de acordo com o modelo constitutivo escolhido pelo usuário. Para resolver a equação do tensor polimério (2.7) usaremos o tensor de conformação A, ^{15), (16), (2}. Este tensor é simétrico e positivo definido, sendo estas importantes propriedades matemática para a construção de matrizes de tranformação ou decomposições. Em geral, a equação para A * é escrita como

$\frac{\partial {\mathbf{A}}^{*}}{\partial t^{*}}+\left({\mathbf{u}}^{*}·{\nabla }^{*}\right){\mathbf{A}}^{*}-\left[{\mathbf{A}}^{\mathbf{*}}{\nabla }^{*}{\mathbf{u}}^{*}+{\nabla }^{*}{{\mathbf{u}}^{\mathbf{*}}}^T{\mathbf{A}}^{*}\right]=\frac{1}{De}{\mathcal{M}}^{*}\left({\mathbf{A}}^{*}\right),$ (2.8)

onde ℳ(A) é definida de acordo com o modelo viscoelástico. As variações nas propriedades reológicas do fluido, como a tensão e viscosidade, são importantes no que diz respeito à escoamentos viscoelásticos. Existem alguns modelos constitutivos que em geral relacionam essas propriedades com o movimento, os quais permitem o estudo do escoamento desses tipos de fluidos. Assim, um modelo constitutivo deve ser usado para descrever o tensor das tensões no escoamento. A tarefa de identificar uma equação reológica adequada para o estudo do escoamento em questão pode não ser tão simples. O comportamento do fluido deve ser expresso por uma equação constitutiva para a tensão T. Uma das maiores dificuldades no estudo da mecânica dos fluidos não Newtonianos está na escolha de um modelo que possa descrever corretamente essas propriedades reológicas. Existem vários modelos constitutivos para este estudo. Neste trabalho, utilizamos o modelo linear de Phan-Thien-Tanner (PTT) ³⁴ para resolver o problema de escoamento eletroosmótico em um microcanal. O modelo PTT é descrito por:

$f\left(tr\right(\mathbf{T}\left)\right)\mathbf{T}+\lambda \overline{\mathbf{T}}=2{\eta }_p\mathbf{D},$ (2.9)

onde $\overline{\mathbf{T}}$ e a função f são dados por:

$\overline{\mathbf{T}}=\frac{\partial \mathbf{T}}{\partial t}+\left(\mathbf{u}·\nabla \right)\mathbf{T}-\left[{\left(\nabla \mathbf{u}\right)}^T·\mathbf{T}+\mathbf{T}·\nabla \mathbf{u}\right]+\xi \left(\mathbf{T}·\mathbf{D}+\mathbf{D}·\mathbf{T}\right),$ (2.10)

$f\left(tr\right(\mathbf{T}\left)\right)=1+\frac{\epsilon \lambda }{{\eta }_p}tr\left(\mathbf{T}\right).$ (2.11)

Portanto, para o modelo PTT, o lado direito da equação (2.7) pode ser reescrito como:

${\mathbf{M}}^{\mathbf{*}}\left({\mathbf{T}}^{\mathbf{*}}\right)=\frac{2(1-\beta )}{Re}{\mathbf{D}}^{\mathbf{*}}-\left(1+\frac{\epsilon ReDe}{1-\beta }\mathrm{tr}\left({\mathbf{T}}^{\mathbf{*}}\right)\right){\mathbf{T}}^{\mathbf{*}}-\xi De\left({\mathbf{T}}^{\mathbf{*}}·{\mathbf{D}}^{\mathbf{*}}+{\mathbf{D}}^{\mathbf{*}}·{\mathbf{T}}^{\mathbf{*}}\right).$ (2.12)

No modelo PTT, a tensão depende do traço de T, e do parâmetro adimensional ε, o qual é relacionado com a viscosidade elongacional para escoamentos desenvolvidos. Outro parâmetro que aparece neste modelo, ξ, essencialmente diz respeito ao movimento não afim entre as moléculas de polímero e o meio contínuo. As cadeias de polímero, dentro do meio, podem se movimentar relativamente à deformação do meio macroscópico, assim cada porção microscópica do polímero pode transmitir apenas uma fração da sua tensão para o meio contínuo ao seu redor. Se $\xi =0$, não existe deslizamento significando que o movimento é afim. Desta maneira, ξ é o parâmetro responsável pela segunda diferença de tensão de cisalhamento ser não nula, levando a fluxos secundários quando o escoamento ocorre em dutos com seção transversal não circular ²⁵. O lado direito da equação do tensor de conformação é dado por:

${\mathcal{M}}^{\mathcal{*}}\left({\mathbf{A}}^{\mathbf{*}}\right)=\left(1+\frac{\epsilon ReDe}{1-\beta }\mathrm{tr}\left({\mathbf{S}}^{\mathbf{*}}\right)\right)\left(\mathbf{I}-{\mathbf{A}}^{\mathbf{*}}\right)-2\xi De\left({\mathbf{B}}^{\mathbf{*}}-{\mathbf{B}}^{\mathbf{*}}{\mathbf{A}}^{\mathbf{*}}\right).$ (2.13)

onde B * é resultante da decomposição do gradiente da velocidade conforme proposto por Fattal and Kupferman ¹⁵. Neste sentido, as equações de movimento são resolvidas para o escoamento regido pelo modelo PTT.

Um grande desafio para pesquisadores da área de computação, especificamente tratando-se de reologia, é resolver a equação (2.7) ou (2.8) para altas viscoelasticidades, ou seja, quanto maior o número de Deborah, maior a dificuldade de o método numérico obter uma solução estável. Com o intuito de suprimir essa instabilidade numérica, Fattal and Kupferman ¹⁵ proporam a reformulação da equação diferencial constitutiva para o logarítmo do tensor de conformação. Numa extensão dessas idéias propostas por Fattal and Kupferman ^{15), (16}, Afonso et al. ² apresentaram uma transformação do tensor de forma genérica, denominada kernel-conformation, que permite aplicar diferentes funções kernel para a matriz de tranformação, na qual a evolução da equação para 𝕜*(A *), pode ser expressada na formulação tensorial como:

$\frac{\partial {\mathbf{𝕜}}^{\mathit{*}}\left({\mathbf{A}}^{\mathbf{*}}\right)}{\partial t^{*}}+\left({\mathbf{u}}^{\mathbf{*}}·{\nabla }^{*}\right){\mathbf{𝕜}}^{\mathit{*}}\left({\mathbf{A}}^{\mathbf{*}}\right)={\mathit{\Omega }}^{\mathit{*}}{\mathbf{𝕜}}^{\mathit{*}}\left({\mathbf{A}}^{\mathbf{*}}\right)-{\mathbf{𝕜}}^{\mathit{*}}\left({\mathbf{A}}^{\mathbf{*}}\right){\mathit{\Omega }}^{\mathit{*}}+2{\mathbb{B}}^{\mathbb{*}}+\frac{1}{De}{\mathbb{M}}^{\mathbb{*}},$ (2.14)

onde 𝔹 and 𝕄 são tensores simétricos construidos pela ortogonalização dos tensores diagonais enquanto que Ω* é um tensor antissimétrico. Assim, a solução numérica pode ser feita através da equação (2.14) ao invés de usar a expressão (2.7). Se consideramos o escoamento newtoniano, o tensor S * na equação de movimento (2.5) é nulo e apenas a velocidade e a pressão são atualizados a cada passo no tempo. Para completar as equações governantes, fica faltando o termo fonte elétrico, que força o fluido a escoar conforme as propriedades elétricas impostas.

Considere um canal plano conforme mostra a figura 1, e que inicialmente o fluido está em repouso dentro do canal. Cargas positivas se acumulam próximo as paredes negativamente carregadas, como mostra a figura 2 e isso faz com que surja um potencial induzido ψ, dando origem a uma fina camada eletricamente carregada, chamada camada de Debye. Em $t=0$ um gradiente de potencial elétrico ∇ϕ é aplicado entre as extremidades de entrada e saida de fluido no canal. Assim, ocorre a movimentação dos íons dentro da fina camada perto das paredes do canal, fazendo com que as partículas neutras sejam arrastadas pela força viscosa, colocando o fluido em movimento. Essencialmente esse é o fenômeno de eletroosmose. A figura 2 mostra o perfil de velocidade eletroosmótico totalmente desenvolvido.

Figura 1:
Ilustração do escoamento eletroosmótico entre duas placas planas paralelas. O fluxo de fluido ocorre entre as duas placas de meia altura igual a H. O comprimento L e largura w são considerados muito maiores que H, mais especificamente para as simulações computacionais $L\approx 20H$.

Figura 2:
As cargas que inicialmente estavam distribuidas de forma homogênea no fluido, se redistribuem dentro do canal, havendo um acúmulo de íons positivamente carregados próximo as paredes negativamente carregadas. Após a aplicação do potencial externo, essas cargas começam a se movimentar próximo as paredes, arrastando as partículas neutras do fluido. O perfil de velocidade mostrado está totalmente desenvolvido para o escoamento eletroosmótico, considerando o gradiente de pressão nulo.

Com o objetivo de resolver a equação do movimento (2.2), temos que incluir o termo forçante elétrico que é dado pela força elétrica

$\mathbf{F}={\rho }_e\mathbf{E},$ (2.15)

onde E é campo elétrico e ρ _e a densidade de carga elétrica. A descrição matemática do fenômeno eletrocinético surge através das equações de Maxwell do eletromagnetismo. Na ausência de campos magnéticos, as equações de Maxwell se tornam:

$\nabla \times \mathbf{E}=0,$ (2.16)

${ϵ}_e\nabla ·\mathbf{E}={\rho }_e,$ (2.17)

onde ε _e é a permissividade elétrica do meio em questão.

Da equação 2.16, podemos escrever $\nabla \times \mathbf{E}=\nabla \times (-\nabla \Phi )=0$, e portanto, usando as expressões 2.16 e 2.17 temos:

$\mathbf{E}=-\nabla \Phi ,$ (2.18)

$\nabla ·{ϵ}_e\nabla \Phi =-{\rho }_e,$ (2.19)

onde o potencial elétrico total Φ e a densidade de carga elétrica ρ _e são funções escalares que em geral dependem das coordenadas espaciais. Note que o potencial Φ é a contribuição total devido as cargas no fluido distribuidas próximo as paredes do canal, o qual chamamos de ψ, e devido ao potencial externo aplicado, que denotaremos por ϕ, ou seja, $\Phi =\varphi +\psi$. Assim, resolvendo a equação (2.19) para o potencial Φ, o campo elétrico E é determinado pela expressão (2.18), e pode ser usado juntamente com a densidade de carga ρ _e para obter o termo fonte da equação de movimento. Portanto além do campo elétrico, devemos determinar ρ _e que é relacionada com a concentração iônica de acordo com a expressão:

${\rho }_e=ez(n^{+}-n^{-}),$ (2.20)

onde e é a constante de carga elementar, $z=z^{+}=z^{-}$ é o número de valência da solução eletrolítica considerando um eletrólito simétrico, e n ⁺ e n ⁻ são as concentrações de íons positivos e negativos respectivamente. Neste sentido, para determinar a densidade de carga ρ ^e , é necessário saber as concentrações n ⁺ e n ⁻ , que são obtidas através da equação de Nernst-Planck, que é dada por:

$\frac{\partial n^{\pm }}{\partial t}=\nabla ·\left(D^{\pm }\nabla n^{\pm }-\mathbf{u}{n}^{\pm }\pm D^{\pm }{n}^{\pm }\frac{ez}{k_{B}T}\nabla \Phi \right),$ (2.21)

onde k _B é a constante de Boltzmann, T a temperatura absoluta, D ⁺ e D ⁻ são os coeficientes de difusão dos íons positivos e negativos respectivamente. Considerando por enquanto um fluido com permissividade εe uniforme e coeficientes de difusão $D^{+}=D^{-}=D$ constantes e, reescrevendo o conjunto de equações (2.19)-(2.21) na forma adimensional, temos:

${\nabla }^{*2}{\Phi }^{*}=-{\rho }_e^{*}$ (2.22)

${\rho }_e^{*}=\delta (n^{*+}-n^{*-}),$ (2.23)

$\frac{\partial n^{*\pm }}{\partial t^{*}}=-{\mathbf{u}}^{*}·{\nabla }^{*}{n}^{*\pm }+\frac{1}{Pe}{\nabla }^{*2}{n}^{*\pm }\pm \frac{1}{Pe}\alpha {\nabla }^{*}·\left(n^{*\pm }{\nabla }^{*}{\Phi }^{*}\right),$ (2.24)

onde $n=n_0{n}^{*}$ e n ₀ a concentração de referência da solução, $\psi ={\zeta }_0{\psi }^{*}$ com ζ ₀ o potencial na parede do canal, $\delta =n_{0}ez{H}^2/{ϵ}_e{\zeta }_0$ um parâmetro adimensional para a simulação, $Pe=UH/D$ o número de Peclet e $\alpha =ez{\zeta }_0/k_{B}T$ é a razão entre a energia potencial e a energia térmica. O conjunto de equações de Poisson-Nernst-Planck adimensionalizadas (2.22) (2.24), que vamos denotar por PNP, são resolvidas numéricamente para obtenção do termo fonte ${\mathbf{F}}^{\mathbf{*}}=-{\rho }_e^{*}{\nabla }^{*}{\Phi }^{*}$ da equação de movimento (2.5).

Podemos simplificar o tratamento das equações PNP, (2.22) (2.24), com o objetivo de obter expressões analíticas para as propriedades do escoamento eletroosmótico n, ψ e ρ _e . Em geral simplificações das equações levam a perda do transiente e é o que ocorre com a equação de transporte iônico (2.24) quando consideramos que a concentração é dada pela distribuição de Boltzmann ¹⁸, ou seja, a solução estacionária da equação (2.24) é dada por:

$n^{*\pm }=\exp (\mp \alpha {\psi }^{*}),$ (2.25)

e substituindo n* ^± na equação (2.23), temos

${\rho }_e^{*}=-2\delta \sinh \left(\alpha {\psi }^{*}\right),$ (2.26)

e para pequenos valores de α, ou seja, energia térmica muito maior que a energia potencial, a equação (2.26) pode ser linearizada, resultando na chamada aproximação de Debye-Hückel (DH):

${\rho }_e^{*}=-2\alpha \delta {\psi }^{*},$ (2.27)

onde é definido o parâmetro adimensional de Debye $\kappa =\sqrt{2\alpha \delta }=H/{\lambda }_D$, com λ _D a largura da camada de Debye.

A descrição do fenômeno através das expressões matemáticas possibilita resolver as equações de movimento para um fluido viscoelástico com um termo forçante que é obitido pela solução das equações PNP. Na próxima seção será descrito o procedimento numérico utilizado.

3 PROCEDIMENTO NUMÉRICO

Para a discretização espacial, foi usado o método de diferenças finitas generalizadas de segunda ordem em uma malha hierárquica, na qual as propriedades como pressão, tensor viscoelástico, potenciais e concentrações elétricas são armazenadas no centro das células, enquanto que as velocidades ficam no centro das faces. As interpolações das propriedades do escoamento são feitas através do método dos mínimos quadrados móveis. Para uma abordagem mais detalhada sobre o sistema HiG-Flow, favor consultar ³¹.

Para o termo difusivo foi usada a discretização de segunda ordem por diferenças centrais. Para a discretização do termo convectivo é usado o método CUBISTA. Para uma abordagem mais detalhada do método, favor consultar ³. Este método de discretização também é usado para os termos convectivos existentes nas equações do tensor viscoelástico e da concentração iônica. A equação de poisson aparece naturalmente quando tentamos resolver o acoplammento pressãovelocidade na equação de movimento usando o método da projeção, seja este incremental ou não incremental. Em todo passo no tempo este tipo de equação precisa ser resolvida. Além disso, a equação de poisson também deve ser resolvida para os potenciais elétricos e também para a concentração iônica nas equações de Poisson-Nernst-Planck. Ao todo são 4 equações de poisson a serem resolvidas. O método da projeção não-incremental ²⁶ foi utilizado para resolver o acoplamento pressão-velocidade. Apesar de a HiG-Flow possibilitar o uso do método incremental, este método de projeção fora escolhido para evitar que oscilações no transiente ocorram durante as simulações ³². Isto não significa problemas quanto a demora na marcha no tempo, pois a própria solução das equações PNP exigem que menores ∆t sejam adotados. Além disso o método de Euler de primeira ordem com o termo difusivo implícito foi utilizado para a integração temporal. Denotaremos o passo no tempo pelo superescrito s, para não confundir com a concentração iônica já definida anteriormente como n. Considere a equação do movimento adimensionalizada (2.5) e, para simplificar a notação, a partir daqui omitiremos o superescrito asterisco, pois só vamos tratar as equações adimensionalizadas. A discretização dessa equação, separando-a nos dois passos devido ao método da projeção é .

$\frac{{\mathbf{u}}^{s+\tau }-{\mathbf{u}}^s}{\Delta t}=-{\mathbf{u}}^s·\nabla {\mathbf{u}}^s+\frac{1}{Re}{\nabla }^2{\mathbf{u}}^{s+\tau }+\nabla ·{\mathbf{S}}^s+{\mathbf{F}}^s,$ (3.1)

$\frac{{\mathbf{u}}^{s+1}-{\mathbf{u}}^{s+\tau }}{\Delta t}=-\nabla p^{s+1},$ (3.2)

onde $s+\tau$ é um tempo intermediário entre os passos s e $s+1$. Aplicando o operador ∇ na equação (3.2) e levando em conta a equação da continuidade,

${\nabla }^2{p}^{s+1}=\frac{1}{\Delta t}\nabla ·u^{s+\tau },$ (3.3)

Assim, a quantidade $p^{s+1}$ é calculada pela solução da equação de poisson e, consequentemente a velocidade no passo $s+1$ é determinada pela substituição de $p^{s+1}$ na equação (3.2). Essencialmente este é o método da projeção não-incremental utilizado. O grande trabalho a ser feito é calcular todos os termos envolvidos na equação (3.1), ou seja, ainda temos que calcular o termo convectivo $-{\mathbf{u}}^s·\nabla {\mathbf{u}}^s$, o termo do tensor viscoelástico $\nabla ·{\mathbf{S}}^s$ e o termo forçante elétrico ${\mathbf{F}}^s={\rho }_e^s\nabla {\Phi }^s$. Como termo difusivo $\frac{1}{Re}{\nabla }^2{\mathbf{u}}^{s+\tau }$ é discretizado no passo $s+\tau$, podemos reescrever a equação (3.1) como sendo

$$\begin{gathered} u_{i,j}^{s+\tau }+\frac{\Delta t}{Re}\left[\frac{2u_{i,j}^{s+\tau }-u_{i+1,j}^{s+\tau }-u_{i-1,j}^{s+\tau }}{\Delta x^2}+\frac{2u_{i,j}^{s+\tau }-u_{i,j+1}^{s+\tau }-u_{i,j-1}^{s+\tau }}{\Delta y^2}\right] \\ =u_{i,j}^s+\Delta t.\mathrm{R}({\mathbf{u}}^s,{\mathbf{S}}^s,p^s,{\mathbf{F}}^s). \end{gathered}$$ (3.4)

onde $\mathrm{R}({\mathbf{u}}^s,{\mathbf{S}}^s,p^s,{\mathbf{F}}^s)$ é a soma dos termos de pressão, convectivo, tensão viscoelástica e força elétrica. A equação (3.4) foi escrita para o compenente u do vetor velocidade $\mathbf{u}=(u,v)$ em duas dimensões. Uma expressão análoga para o componente v é resolvida no mesmo ”laço for”que percorre todas as células do domínio. Portanto, determinando a função $\mathrm{R}({\mathbf{u}}^s,{\mathbf{S}}^s,p^s,{\mathbf{F}}^s)$, o sistema de equações lineares pode ser resolvido numéricamente. O termo correspondente a pressão é atualizado de forma que $p^{s+1}=p^s+dp$, com d _p calculado na equação (3.3).

Ao invés de resolver diretamete a equação do tensor viscoelástico T, foi usado o kernelconformation, sendo discretizada a expressão para o tensor 𝕜, equação (2.14). Essencialmente, o sistema HiG-Flow possui uma rotina pra calcular o tensor kernel 𝕜, outra rotina para marchar esse tensor kernel no tempo e por fim fazer a transformação do kernel atualizado para obter o tensor polimérico T e o tensor S que deve ser usado na equação de movimento. A discretização para a marcha no tempo do tensor kernel é feita de forma explícita, ou seja,

${\mathbf{𝕜}}^{s+1}={\mathbf{𝕜}}^s+\Delta t\left[-\left({\mathbf{u}}^s·\nabla \right){\mathbf{𝕜}}^s+{\mathit{\Omega }}^s{\mathbf{𝕜}}^s-{\mathbf{𝕜}}^s{\mathit{\Omega }}^s+2{\mathbb{B}}^s+\frac{1}{De}{\mathbb{M}}^s\right],$ (3.5)

Uma rotina é utilizada para calcular o termo ${\mathit{\Omega }}^s{\mathbf{𝕜}}^s-{\mathbf{𝕜}}^s{\mathit{\Omega }}^s$ e outra para obter $2{\mathbb{B}}^s+\frac{1}{De}{\mathbb{M}}^s$. Resta o termo convectivo, $\left({\mathbf{u}}^s·\nabla \right){\mathbf{𝕜}}^s$ que é calculado usando o método CUBISTA. Este método garante maior estabilidade a solução numérica quando por exemplo comparado ao método upwind de diferenças centrais.

Além da influência das tensões viscoelásticas do modelo PTT, a existência do termo fonte elétrico impõe um perfil de velocidade característico no escoamento. A discretização dos termos eletrocinéticos foi realizada para que o código funcione no HiG-Flow tanto para o fenômeno físico completo descrito pelo conjunto de equações de Poisson-Nernst-Planck (PNP), quanto para a aproximação de Debye-Hückel (DH). Vamos começar com a equação de transporte iônico (2.24). Omitindo o superescrito ±, a equação para as concentrações pode ser reescrita como sendo:

$\frac{n^{s+1}-n^s}{\Delta t}=-{\mathbf{u}}^s·\nabla n^s+\frac{1}{Pe}{\nabla }^2{n}^{s+1}\pm \frac{1}{Pe}\alpha \nabla ·\left(n^s\nabla {\Phi }^s\right),$ (3.6)

onde ${\Phi }^s=\varphi +\psi$. Assim como na equação (3.1), o termo difusivo é calculado implicitamente. Desta forma,

$$\begin{gathered} n_{i,j}^{s+1}+\frac{\Delta t}{Pe}\left(\frac{2n_{i,j}^{s+1}-n_{i-1,j}^{s+1}-n_{i+1,j}^{s+1}}{\Delta x^2}+\frac{2n_{i,j}^{s+1}-n_{i,j-1}^{s+1}-n_{i,j+1}^{s+1}}{\Delta y^2}\right) \\ =n^s+\Delta t\mathrm{N}(n^s,{\mathbf{u}}^s,{\Phi }^s), \end{gathered}$$ (3.7)

onde

$\mathrm{N}(n^s,{\mathbf{u}}^s,{\Phi }^s)=-{\mathbf{u}}^s·\nabla n^s\pm \frac{1}{Pe}\alpha \nabla ·\left(n^s\nabla {\Phi }^s\right).$ (3.8)

O termo convectivo, primeiro à direita da igualdade (3.8), é obtido a cada passo no tempo pelo método CUBISTA. Já o termo potencial é discretizado por diferenças centrais:

$$\begin{gathered} \pm \frac{1}{Pe}\alpha \nabla ·\left(n^s\nabla {\Phi }^s\right)\approx \pm \frac{1}{Pe}\alpha \frac{(n_{i+1,j}^s-n_{i-1,j}^s)({\varphi }_{i+1,j}^s-{\varphi }_{i-1,j}^s+{\psi }_{i+1,j}^s-{\psi }_{i-1,j}^s)}{2\Delta x}\mathrm{\backslash nonumber} \\ \pm \frac{1}{Pe}\alpha \frac{(n_{i,j+1}^s-n_{i,j-1}^s)({\varphi }_{i,j+1}^s-{\varphi }_{i,j-1}^s+{\psi }_{i,j+1}^s-{\psi }_{i,j-1}^s)}{2\Delta y} \\ \pm \frac{1}{Pe}\alpha \delta n_{i,j}^s(n_{i,j}^{s+}-n_{i,j}^{s-}). \end{gathered}$$ (3.9)

Assim, n ^± são determinados pela substituição de $\mathrm{N}(n^s,{\mathbf{u}}^s,{\Phi }^s)$ em na equação (3.7). Em seguida, o potencial ψ precisa ser atualizado pela equação de poisson

$\frac{{\psi }_{i-1,j}^{s+1}-2{\psi }_{i,j}^{s+1}+{\psi }_{i+1,j}^{s+1}}{\Delta x^2}+\frac{{\psi }_{i,j-1}^{s+1}-2{\psi }_{i,j}^{s+1}+{\psi }_{i,j+1}^{s+1}}{\Delta y^2}=-\delta (n^{(s+1)+}-n^{(s+1)-}),$ (3.10)

a cada passo de tempo, visto que ele está diretamente acoplado com a concentração iônica. Com o valor de n ^± e o valor de ∇Φ calculados, o termo fonte elétrico é determinado por $\mathbf{F}=-{\rho }_e\nabla \Phi$, que pode ser escrito como:

$\mathbf{F}=-\delta (n_{i,j}^{s+}-n_{i,j}^{s-})\left[\begin{array}{c}\frac{({\varphi }_{i+1,j}^s-{\varphi }_{i-1,j}^s+{\psi }_{i+1,j}^s-{\psi }_{i-1,j}^s)}{2\Delta x}\\ \frac{({\varphi }_{i,j+1}^s-{\varphi }_{i,j-1}^s+{\psi }_{i,j+1}^s-{\psi }_{i,j-1}^s)}{2\Delta y}\end{array}\right]$ (3.11)

3.1 Condições de contorno nas paredes do canal

Nas paredes do canal, condições de Neumann são aplicadas para a pressão p e para o potencial eletrostático ϕ enquanto que ψ assume o valor potencial de referência ζ0. As concentrações nas paredes são dadas pela distribuição de Boltzmann (3.15) e para a velocidade é imposta a condição de não deslizamento.

$\hat{\mathbf{n}}·\nabla p=0,$ (3.12)

$\hat{\mathbf{n}}·\nabla \varphi =0,$ (3.13)

$\psi ={\zeta }_0,$ (3.14)

$n^{\pm }=n_0\exp (\mp \alpha {\zeta }_0),$ (3.15)

$\mathbf{u}=\mathbf{0},$ (3.16)

3.2 Condições de contorno na entrada e saída do canal

Para a pressão, condições de Dirichlet são aplicadas, assim a entrada e saída do canal tem exatamente a mesma pressão, resultando em um escoamento eletrosmótico puro, isto é, com p _in igual a p _out . Para o potencial eletrostático ϕ também impomos os valores na entrada e saída do canal, com ${\varphi }_{in}\ne {\varphi }_{out}$ obtendo assim um $\nabla \varphi \ne 0$. A concentração assume o valor de referência n ₀ na entrada do canal enquanto que na saída a condição de Neumann $\hat{\mathbf{n}}·\nabla n^{\pm }=0$ é imposta para que o escoamento seja totalmente desenvolvido. Impomos condição de Neumann para a velocidade tanto na entrada como na saída do canal, para que o escoamento seja desenvolvido apenas sob a ação da força elétrica.

Abaixo é descrita a sequência correta de cálculos para obter a solução em cada passo de tempo s. Para $s=0$, temos:

• passo 1: impor as condições de contorno adequadas para a geometria utilizada.

• passo 2: calcular a velocidade intermediária u ^τ usando a equação (3.4). Note que todas as propriedades, p ⁰, u ⁰, S ⁰, ψ ⁰, ϕ ⁰, n ^0(±) devem ser inicializadas no volume, de acordo com o condição inicial desejada. Por exemplo, para o escoamento entre placas paralelas $p^0=0,{\mathbf{u}}^0=\mathbf{0},{\mathbf{S}}^0=\mathbf{0},{\psi }^0=0,{\varphi }^0=0,n^{0(\pm )}=1$ .

• passo 3: calcular o termo de pressão, p ¹ resolvendo a equação de poisson (3.3).

• passo 4: calcular a velocidade atualizada u ¹ usando a equação (3.2).

• passo 5: calcular 𝕜¹ usando a equação (3.5) e realizar a transformação para obter ∇ · S ¹.

• passo 6: resolver a equação de transporte iônica (3.7) para obter n ^1(±) .

• passo 7: calcular o potencial ψ ¹ resolvendo a equação de poisson (3.10).

• passo 8: calcular o potencial ϕ ¹ resolvendo a equação de poisson ${\nabla }^2\varphi =0$ .

• passo 9: calcular o termo fonte elétrico F ¹, dado por (3.11).

O procedimento acima deve ser repetido para cada passo de tempo até atingir a solução estacionária, se existir. Para o caso da aproximação de Debye-Hückel, o passo 6 não é executado.

4 RESULTADOS

Um dos pontos que vale ressaltar quando queremos simular escoamentos eletroosmóticos é o refinamento da malha perto das paredes do canal. Se a malha for muito grossa naquela região, informações podem ser perdidas e o resultado final deve ser afetado. O efeito do refinamento da malha é mostrado na figura 3a. Note que para uma distância mínima entre dois pontos vizinhos de $2,5\times {10}^{-2}$ existe um grande erro em comparação com a curva analítica. Para a malha um pouco mais refinada, $\Delta y/H=3,125\times {10}^{-3}$, a curva se aproxima do resultado analítico, assim como para a malha mais refinada com $\Delta y/H=7,8215\times {10}^{-4}$. Para as simulações com valores considerados altos para o parâmetro de Debye, $\kappa >100$ esse nível de refinamento é necessário para que se obtenha um número razoável de pontos próximo as paredes. A figura 3b ilustra a malha utilizada nas simulções que exigem grande refinamento perto das paredes, onde ocorrem variações significativas das propriedades do escoamento. Esta malha possui 7 níveis de refinamento com $\Delta y_{min}/H=7,8215\times {10}^{-4}$.

Figura 3:
a) Efeito do refinamento próximo as paredes do canal e b) malha utilizada nas simulações de escoamento entre placas paralelas.

Nas paredes do canal, condições de Neumann são aplicadas para a pressão p e para o potencial eletrostátic ϕ enquanto que ψ assume o valor potencial de referência ζ ₀. As concentrações nas paredes são dadas pela distribuição de Boltzmann, $n^{\pm }=n_0\exp (\mp \alpha {\zeta }_0)$ e, para a velocidade é imposta a condição de não deslizamento. Para a pressão, condições de Dirichlet são aplicadas, assim a entrada e saída do canal tem exatamente a mesma pressão, resultando em um escoamento eletroosmótico puro, isto é, com p _in igual a p _out . Para o potencial eletrostático ϕ também impomos os valores na entrada e saída do canal, com ${\varphi }_{in}\ne {\varphi }_{out}$ obtendo assim um $\nabla \varphi \ne 0$. A concentração assume o valor de referência n ₀ na entrada do canal enquanto que na saída a condição de Neumann $\hat{\mathbf{n}}·\nabla n^{\pm }=0$ é imposta para que o escoamento seja totalmente desenvolvido. Impomos condição de Neumann para a velocidade tanto na entrada como na saída do canal, para que o escoamento seja desenvolvido apenas sob a ação da força elétrica.

4.1 Fluido newtoniano

No caso de um fluido newtoniano, o efeito da força elétrica torna o perfil de velocidade dependente da distância até a parede do canal e da espessura da camada de Debye, isto é, do parâmetro κ. Quanto menor a espessura da camada de Debye, maior é a concentração de cargas perto das paredes, onde fortes variações do perfil de velocidade devem acontecer.

A figura 4a mostra o perfil de velocidade para diferentes espessuras da camada de Debye, e a figura 4b amplia a imagem próximo a parede do canal. A medida que diminui a espessura da camada de Debye, ou seja, aumentando o parâmetro κ, ocorrem variações mais acentuadas no perfil de velocidade perto das paredes. As curvas estão normalizadas pela velocidade de Helmholtz-Smoluchowski, $u_{sh}=-{ϵ}_e{\zeta }_0{E}_x/\eta$. Para estas simulações, $Re=0,001,Pe=1,0$. As concentrações iônicas n ^± são mostradas na figura 5a, e na ampliação da região próxima a parede do canal, figura 5b, nota-se perfeitamente a maior concentração de cargas perto da parede a medida que o parâmetro de Debye cresce de $\kappa =10\text{a}300$. O acúmulo de cargas perto da parede faz com que o potencial ψ varie nesta região. Para $\kappa =10$ mudanças no potencial são sentidas a maiores distâncias da parede quando comparada a curva mostrada para $\kappa =20$. A medida que a camada de Debye se torna mais fina, $\kappa =50$, o potencial varia rapidamente perto das paredes do canal e consequentemente a concentração de cargas tende ficar cada vez maior naquela região como mostrado para $\kappa =300$, figuras 6 e 5.

Figura 4:
a) Perfil de velocidade para escoamento newtoniano entre placas paralelas com diferentes valores do parâmetro de Debye κ e em b) ampliação da região próxima a parede do canal.

Figura 5:
a) Concentração iônica n ⁺ e n ⁻ para escoamento newtoniano entre placas paralelas com diferentes valores do parâmetro de Debye κ e em b) ampliação da região próxima a parede do canal.

Figura 6:
a) Potencial eletroosmótico ψ para escoamento newtoniano entre placas paralelas com diferentes valores do parâmetro de Debye κ e em b) ampliação da região próxima a parede do canal.

Para o fluido viscoelástico, o modelo PTT foi usado para resolver este problema. Primeiramente apresentamos os resultados para o model PTT completo e depois para o modelo simplificado sPTT.

4.2 Modelo de Phan-Thien/Tanner-PTT

O perfil de velocidade depende da distância à parede, y, do parâmetro de Debye κ e também do tempo de relaxação do fluido que é relacionado com o número de Deborah $D{e}_{\kappa }=\lambda \kappa u_{sh}$, como mostrado na figura 7. As simulações foram realizadas para $Re=0,001$ e $Pe=1,0$.

Figura 7:
Perfil de velocidade para escoamento viscoelástico entre placas paralelas para $\kappa =10,\epsilon =\xi =0,01$.

Para $D{e}_{\kappa }=0,5$, o comportamento do fluxo se aproxima ao de um escoamento newtoniano. Se $D{e}_{\kappa }=2,0$ nota-se que o perfil de velocidade sofre maiores alterações em relação ao perfil newtoniano, e isso ocorre devido a baixas viscosidades perto das paredes do canal, influenciados pelo aparecimento de tensões cisalhantes do modelo PTT. As tensões normais e cisalhantes são mostradas na figura 8a e 8b respectivamente. Os resultados são similares aos encontrados na literatura. Para detalhes sobre a solução analítica para o estado estacionário, favor consultar ¹².

Figura 8:
Tensões para o fluido PTT em função do número de Deborah, para $\kappa =10$ e $\kappa =100,\epsilon =\xi =0,01$. a) Tensão T _xx e em b) T _xy .

4.3 Modelo de Phan-Thien/Tanner simplificado-sPTT

Nesse modelo $\xi =0$ na expressão (2.13). Além disso, ao invés de considerarmos o escoamento eletroosmótico puro, vamos impor um gradiente de pressão não nulo no canal. Neste sentido, as soluções para a velocidade e o tensor viscoelástico no estado estacionário para o modelo sPTT dependem do parâmetro Γ definido como:

$\Gamma =\frac{H^2}{{ϵ}_e{\zeta }_0{E}_x}\frac{dp}{dx}.$ (4.1)

Neste sentido, para um gradiente de pressão maior que zero, $dp/dx>0$, o fluido é ”empurrado”no sentido contrário ao escoamento conforme mostra a figura 9 para $\mathrm{\Gamma }=1,0,\mathrm{\Gamma }=2,0\mathrm{e}\mathrm{\Gamma }=2,778$, representados triângulos apontando para baixo, losangos e triângulos apontando para esquerda respectivamente. Por outro lado, se dp/dx < 0, o fluido é empurrado no sentido do escoamento, resultando nos perfis mostrados na figura 9 para $\mathrm{\Gamma }=-0,5\mathrm{e}\mathrm{\Gamma }=-1,778$, representados por círculos e quadrados respectivamente. Para detalhes sobre a solução analítica para o estado estacionário deste problema, favor consultar ¹.

Figura 9:
Perfil de velocidade para escoamento newtoniano entre placas paralelas para $\kappa =20$. As curvas mostram o efeito do gradiente de pressão aplicado, relacionado com o parâmetro Γ.

A figura 10a mostra o perfil de velocidade para o escoamento eletroosmótico do fluido sPTT para $\mathrm{\Gamma }=-1,0$ com triângulos apontando para direita e $\mathrm{\Gamma }=-2,77$ com triângulos apontando para a esquerda. Nota-se que assim como no caso do fluido PTT, o perfil de velocidade no modelo sPTT é influenciado pelo número de Deborah, que nestas simulações é $D{e}_{\kappa }=2,5$. Este efeito é melhor visualizado para as curvas mostradas quando $\mathrm{\Gamma }=0$ na figura 10a, representadas pelas linhas sólidas para $D{e}_{\kappa }=0$ e para $D{e}_{\kappa }=2,5$, onde nota-se claramente a mudança relativa ao perfil newtoniano. As tensões normais e de cisalhamento sofrem um deslocamento em relação àquelas que correspondem ao escoamento eletroosmótico puro, no qual $\mathrm{\Gamma }=0$ representado pelas linhas sólidas, como indica a figura 10b. Para o componente normal T _xx , as curvas indicam uma diminuição da tensão perto da parede do canal quando o gradiente de pressão aplicado é positivo, $\mathrm{\Gamma }=2,77$, representado por triângulos apontando para cima. Por outro lado, se o gradiente de pressão for negativo, $\mathrm{\Gamma }=-1,0$, então T _xx aumenta discretamente perto das paredes do canal, como mostra a curva representada por triângulos apontando para baixo. Comportamento similar é observado para a magnitude da tensão de cisalhamento T _xy .

Figura 10:
Efeito do gradiente de pressão aplicado no escoamento eletroosmótico de fluido viscoelástico entre placas paralelas. Simulações realizadas com $\kappa =20,D{e}_{\kappa }=2,5\text{e}\mathrm{\Gamma }=-1,\mathrm{\Gamma }=2,77$. Em a) velocidade e em b) componentes do tensor polimérico.

5 CONCLUSÃO

Foram estudados escoamentos eletroosmóticos entre placas planas. Foram utilizados os modelos de Poisson-Nernst-Planck e de Debye-Hückel. No que diz respeito as simulações nos canais de placas planas paralelas, este foi o início da resolução dos problemas eletroosmóticos que podem ser feitas com a HiG-Flow. Como o código ainda está sendo desenvolvido, fizemos todas as verificações dos resultados obtidos. Foi estudado o efeito do campo elétrico aplicado em um escoamento regido pelo modelo constituvo de PTT e realizadas comparações com resultados analíticos e numéricos já existentes na literatura para o problema. Malhas adaptativas foram usadas devido ao refinamento necessário perto das paredes eletricamente carregadas. Este é um dos principais problemas em relação ao tempo de simulação para escoamentos sujeitos a forças elétricas. Na verdade, esta é a solução para o problema, pois usar uma malha regular seria totalmente inviável tanto em termos de memória quanto em tempo de máquina, principalmente se desejar usar $\kappa \ge 100$. Outro problema que ocorre nas simulações é quando se tem altos números de Deborah. Para cada tipo de modelo viscoelástico, geometria de escoamento existe um número de Deborah crítico para o qual as simulações podem ser fortemente perturbardas. Para os escoamentos estudados aqui, instabilidades no método tornam-se mais evidentes para $D{e}_{\kappa }>2$. Um estabilador aplicado por meio do ”kernel-conformation” ajuda a melhorar a estabilidade da simulação, mas quanto maior for De maiores perturbações podem deixar o resultado insatisfatório.

AGRADECIMENTOS

Os autores agradecem o apoio financeiro dado pela Fapesp pelo projeto CEPID processo 2013/07375-0. A. Castelo agradece o processo CNPq 307483/2017-7. W. S. Bezerra agradece o suporte dado pela FACET-UFGD (Faculdade de Ciências Exatas e Tecnologia da Universidade Federal da Grande Dourados).

REFERÊNCIAS

• ¹
  A. Afonso, M. Alves & F. Pinho. Analytical solution of mixed electro-osmotic/pressure driven flows of viscoelastic fluids in microchannels. Journal of Non-Newtonian Fluid Mechanics, 159(1-3) (2009), 50-63.
• ²
  A. Afonso, F. Pinho &M. Alves . The kernel-conformation constitutive laws. Journal of Non-Newtonian Fluid Mechanics , 167 (2012), 30-37.
• ³
  M. Alves , P. Oliveira & F. Pinho. A convergent and universally bounded interpolation scheme for the treatment of advection. International journal for numerical methods in fluids, 41(1) (2003), 47-75.
• ⁴
  S. Balay, S. Abhyankar, M.F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W.D. Gropp, D. Kaushik, M.G. Knepley, L.C. McInnes, K. Rupp, B.F. Smith, S. Zampini, H. Zhang & H. Zhang . PETSc Web page (2017). URL http://www.mcs.anl.gov/petsc
  » http://www.mcs.anl.gov/petsc
• ⁵
  F. Bianchi, R. Ferrigno & H. Girault. Finite element simulation of an electroosmotic-driven flow division at a T-junction of microscale dimensions. Analytical Chemistry, 72(9) (2000), 1987-1993.
• ⁶
  H. Bruus. “Theoretical microfluidics. Oxford master series in condensed matter physics”. Oxford University Press, Oxford, UK (2008).
• ⁷
  D. Burgreen & F. Nakache. Electrokinetic flow in ultrafine capillary slits1. The Journal of Physical Chemistry, 68(5) (1964), 1084-1091.
• ⁸
  A. Castelo, A. Afonso & W. Souza. A finite difference method in hierarquical grids for viscoelastic fluid flow simulations. (2018). In preparation.
• ⁹
  M.L. Chabinyc, D.T. Chiu, J.C. McDonald, A.D. Stroock, J.F. Christian, A.M. Karger & G.M. Whitesides. An integrated fluorescence detection system in poly (dimethylsiloxane) for microfluidic applications. Analytical Chemistry , 73(18) (2001), 4491-4498.
• ¹⁰
  D.L. Chapman. LI. A contribution to the theory of electrocapillarity. The London, Edinburgh, and Dublin philosophical magazine and journal of science, 25(148) (1913), 475-481.
• ¹¹
  P. Debye & E. Hückel. De la theorie des electrolytes. I. abaissement du point de congelation et phenomenes associes. Physikalische Zeitschrift, 24(9) (1923), 185-206.
• ¹²
  S. Dhinakaran, A. Afonso , M. Alves & F. Pinho. Steady viscoelastic fluid flow between parallel plates under electro-osmotic forces: Phan-Thien-Tanner model. Journal of colloid and interface science, 344(2) (2010), 513-520.
• ¹³
  E.A. Doherty, R.J. Meagher, M.N. Albarghouthi & A.E. Barron. Microchannel wall coatings for protein separations by capillary and chip electrophoresis. Electrophoresis, 24(1-2) (2003), 34-54.
• ¹⁴
  P. Dutta & A. Beskok. Analytical solution of combined electroosmotic/pressure driven flows in twodimensional straight channels: finite Debye layer effects. Analytical chemistry, 73(9) (2001), 1979- 1986.
• ¹⁵
  R. Fattal & R. Kupferman. Constitutive laws for the matrix-logarithm of the conformation tensor. Journal of Non-Newtonian Fluid Mechanics , 123(2) (2004), 281-285.
• ¹⁶
  R. Fattal &R. Kupferman . Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation. Journal of Non-Newtonian Fluid Mechanics , 126(1) (2005), 23-37.
• ¹⁷
  R.A. Finkel & J.L. Bentley. Quad trees a data structure for retrieval on composite keys. Acta informatica, 4(1) (1974), 1-9.
• ¹⁸
  M. Fixman. The Poisson-Boltzmann equation and its application to polyelectrolytes. The Journal of Chemical Physics, 70(11) (1979), 4995-5005.
• ¹⁹
  M. Gouy. Sur la constitution de la charge électrique à la surface d’un électrolyte. J. Phys. Theor. Appl., 9(1) (1910), 457-468.
• ²⁰
  H.V. Helmholtz. Studien über electrische Grenzschichten. Annalen der Physik, 243(7) (1879), 337-382.
• ²¹
  J. Lin, L.M. Fu & R.J. Yang. Numerical simulation of electrokinetic focusing in microfluidic chips. Journal of Micromechanics and Microengineering, 12(6) (2002), 955.
• ²²
  H. Park & W. Lee. Helmholtz-Smoluchowski velocity for viscoelastic electroosmotic flows. Journal of colloid and interface science , 317(2) (2008), 631-636.
• ²³
  N.A. Patankar & H.H. Hu. Numerical simulation of electroosmotic flow. Analytical Chemistry , 70(9) (1998), 1870-1881.
• ²⁴
  R. Peng & D. Li. Effects of ionic concentration gradient on electroosmotic flow mixing in a microchannel. Journal of colloid and interface science , 440 (2015), 126-132.
• ²⁵
  N. Phan-Thien. A nonlinear network viscoelastic model. Journal of Rheology, 22(3) (1978), 259-283.
• ²⁶
  A. Quarteroni, F. Saleri & A. Veneziani. Factorization methods for the numerical approximation of Navier-Stokes equations. Computer methods in applied mechanics and engineering, 188(1-3) (2000), 505-526.
• ²⁷
  F.F. Reuss. Sur un nouvel effet de l’électricité galvanique. Mem. Soc. Imp. Natur. Moscou, 2 (1809), 327-337.
• ²⁸
  C. Rice & R. Whitehead. Electrokinetic flow in a narrow cylindrical capillary. The Journal of Physical Chemistry , 69(11) (1965), 4017-4024.
• ²⁹
  M.v. Smoluchowski. Versuch einer mathematischen Theorie der Koagulationskinetik kolloider Lo¨sungen. Zeitschrift für physikalische Chemie, 92(1) (1918), 129-168.
• ³⁰
  L. Song, L. Yu, Y. Zhou, A.R. Antao, R.A. Prabhakaran & X. Xuan. Electrokinetic instability in microchannel ferrofluid/water co-flows. Scientific reports, 7 (2017), 46510.
• ³¹
  F.S. Sousa, C.F.A. Lages, J.L. Ansoni, A. Castelo & A. Simao. A finite difference method with meshless interpolation for fluid flow simulations in hierarchical grids. Journal of Computational Physics, (2019).
• ³²
  F.S. Sousa , C.M. Oishi & G.C. Buscaglia. Spurious transients of projection methods in microflow simulations. Computer methods in applied mechanics and engineering , 285 (2015), 659-693.
• ³³
  G. Tang, X. Li, Y. He & W. Tao. Electroosmotic flow of non-Newtonian fluid in microchannels. Journal of Non-Newtonian Fluid Mechanics , 157(1-2) (2009), 133-137.
• ³⁴
  N.P. Thien & R.I. Tanner. A new constitutive equation derived from network theory. Journal of NonNewtonian Fluid Mechanics, 2(4) (1977), 353-365.
• ³⁵
  C. Yang &D. Li . Analysis of electrokinetic effects on the liquid flow in rectangular microchannels. Colloids and surfaces A: physicochemical and engineering aspects, 143(2-3) (1998), 339-353.
• ³⁶
  C. Zhao, E. Zholkovskij, J.H. Masliyah &C. Yang . Analysis of electroosmotic flow of power-law fluids in a slit microchannel. Journal of colloid and interface science , 326(2) (2008), 503-510.

Datas de Publicação

• Publicação nesta coleção
  30 Nov 2020
• Data do Fascículo
  Sep-Dec 2020

Histórico

• Recebido
  29 Ago 2018
• Aceito
  06 Ago 2020