UNIVERSIDADE ESTADUAL DE CAMPINAS
FACULDADE DE ENGENHARIA MECÂNICA
DEPARTAMENTO DE MECÂNICA COMPUTACIONAL
TRABALHO DE GRADUAÇÃO EM ENGENHARIA DE CONTROLE E AUTOMAÇÃO
Daniel Scalioni Carvalho
Estudo da simulação de fluidos com o método DPD
Campinas / SP
2013
Daniel Scalioni Carvalho
Estudo da simulação de fluidos com o método DPD
Trabalho de Graduação apresentado como exigência
parcial para obtenção do Diploma de Graduação em
Engenharia de Controle e Automação da Universidade
Estadual de Campinas.
Orientador: Prof. Dr. Luiz Otávio Saraiva Ferreira
Campinas / SP
2013
Agradecimentos
Primeiramente agradeço a Deus por ter me dado saúde, inteligência, paciência e
perseverança não só para finalizar este trabalho de graduação, mas também, por todo o
percurso acadêmico que realizei nesta tão almejada Universidade.
Minha mãe, Natalia Scalioni, e, meu pai, José C. Carvalho, que sempre me instigaram a
ter fome de conhecimento e a não me contentar com pouco, a não ser uma das melhores
universidades do país.
Agradeço minha noiva Débora C. Cantador que sempre esteve ao meu lado, me
apoiando em qualquer que seja a ocasião, tanto na vida da faculdade quanto fora dela. Sem
mencionar sua família que me apoiavam a cada passo mais duro que eu queria dar.
Os amigos da faculdade também foram essenciais, em especial minha cavalaria de
estudos, integrada por Leonardo Reis e Leonardo Junqueira, com algumas ajudas essenciais
de Conrado Miranda.
Sem se esquecer do grupo de pesquisas, liderados pelo professor Luiz Otávio Saraiva
Ferreira, que a qualquer dúvida ou problemas de software estavam ali prontos para
compartilhar suas experiências para que todos saíssem ganhando.
Mas ai vai um agradecimento especial ao professor, orientador e amigo que me
acompanhou por 3 anos, desde o dia em que bati à porta dele, como um novato, querendo
participar de uma iniciação científica.
E finalmente, agradeço a universidade que sempre me forneceu todos os recursos
necessários para a melhor formação que eu poderia ter.
3
Sumário
AGRADECIMENTOS ......................................................................................................... 3
RESUMO ..................................................................................................................... 4
1.
INTRODUÇÃO ....................................................................................................... 5
2.
FUNDAMENTAÇÃO TEÓRICA ..................................................................................... 7
3.
MÉTODOS ......................................................................................................... 27
4.
RESULTADOS E DISCUSSÕES ................................................................................... 32
5.
TRABALHOS FUTUROS .......................................................................................... 41
6.
CONCLUSÃO....................................................................................................... 41
7.
BIBLIOGRAFIA ..................................................................................................... 43
ANEXO A – ARQUIVO DE CONFIGURAÇÃO DO HOOMD-BLUE................................................... 48
ANEXO B – PROGRAMA EM MATLAB ............................................................................... 50
Resumo
Este trabalho trata da validação do simulador Hoomd-blue na simulação de fluidos em
microcanais usando o método DPD (Dissipative Particle Dynamics). A diferença do Hoomdblue é a utilização intensiva de GPGPUs (General-Purpose Graphics Processing Units) na
aceleração do processamento dos cálculos massivos necessários em dinâmica de partículas.
O trabalho consiste em gerar um script de simulação que reproduz o deslocamento de um
fluido num tubo de seção quadrada, e assim, evidenciar, inicialmente, qualitativamente o
perfil parabólico de velocidades. Com isso, será possível obter o perfil de velocidades das
partículas ao longo de uma seção do tubo e compará-lo com resultados analíticos e da
própria literatura.
Palavras-chave: SIMULAÇÃO, DPD, DINÂMICA MOLECULAR, GPGPU, SISTEMAS DE
PARTÍCULAS.
4
1. Introdução
1.1. Apresentação
A miniaturização de dispositivos é um processo extremamente comum no âmbito da
eletrônica. Visa minimizar perdas em espaço físico, potência utilizada e principalmente de
recursos de manufatura, deixando-os mais atrativos para sua comercialização. Com esse
mesmo intuito, outros tipos de dispositivos estão sendo miniaturizados, e são
genericamente chamados de MEMS (Micro-Eletro-Mechanical Systems) [1].
O universo dos dispositivos MEMS é extremamente vasto, com aplicações nas áreas de
telecomunicações, automobilística, medicina, farmacêutica, química, etc. Nas áreas de
química e bioquímica, o interesse maior é em dispositivos microfluidicos tais como reatores
[2] [3] ou controladores de formação de gota [4], dentre outras aplicações [5].
Estes dispositivos são um subconjunto do MEMS chamados Lab-on-a-chip [6]. Esta
nova tecnologia representa uma revolução nos experimentos laboratoriais. Atualmente,
estas aplicações são as que mais usufruem das descobertas em microfluidica. O
desenvolvimento e aperfeiçoamento destes microlaboratórios se baseiam em simulações de
microfluidos. Existem vários métodos de simulação fluidica [7] [8], no entanto, foi nesta
década que os métodos de partículas [9] se tornaram mais acessíveis apresentando
importantes vantagens para microdispositivos.
A utilização de métodos de partículas era restrita devido ao longo tempo de
processamento requerido em computadores convencionais. No entanto, hoje é uma das
áreas que mais há pesquisas em métodos numéricos, pois, com o advento de processadores
massivamente paralelos, o tempo de execução reduziu-se a patamares competitivos com os
métodos convencionais. Estes processadores são as chamadas GPGPU [10] que, utilizando-se
do paradigma de processamento e programação paralelos, conseguem ter desempenhos
muito superiores que das CPUs a custos e consumo de energia muito menores.
Um dos problemas da comunidade científica na utilização de todo esse potencial é a
migração para este novo paradigma. Para resolver isso que o Hoomd-blue está sendo
desenvolvido. No seu estado de desenvolvimento atual o Hoomd-blue roda em uma única
GPU NVIDIA, e é programado na linguagem CUDA C. Seu desempenho é superior a um
5
simulador semelhante rodando em um cluster de CPUs com qualquer número de núcleos
maior que 32, quantidade, a partir da qual, os custos de comunicação pioram o
desempenho. Está também em desenvolvimento uma versão multi-GPU.
1.2. Motivação
O Hoomd-blue é utilizado no Laboratório de Simulações Multifísicas, dirigido pelo meu
orientador, Prof. Dr. Luiz Otávio Ferreira, para a simulação de flúidos em microcanais, com
potencial de aplicação em problemas das áreas médica [11], farmacêutica [5], química [12],
dentre outras.
Esse software está sendo utilizado para que o grupo ganhe experiência na modelagem
de problemas de microfluidica com o método DPD, enquanto um simulador dedicado está
sendo desenvolvido para execução em cluster de GPUs. Portanto é importante que os
resultados obtidos com o Hoomd-blue sejam validados, pois serão utilizados posteriormente
para validação do novo simulador.
1.3. Objetivos
1.3.1. Objetivo Geral
Validar o resultado gerado pelo software Hoomd-blue, utilizando o método DPD, na
simulação de um sistema fluidico no interior de um tubo retangular com base em literaturas
e resultados analíticos.
1.3.2. Objetivos Específicos
1. Configurar um ambiente para a utilização do software Hoomd-blue.
2. Programar o software para realizar a simulação proposta.
3. Avaliar a resposta gerada pelo software do perfil de velocidades do fluido
tomando-se como base soluções analíticas consagradas da literatura bem como
resultados oriundos de outros simuladores padrões na comunidade.
6
2. Fundamentação Teórica
2.1. Microfluidica
Fluido é uma substância que se deforma continuamente sob a aplicação de uma
tensão de cisalhamento (tangencial), não importa quão pequena ela seja [13]. No entanto,
esta definição é muito abrangente para os fluidos tratados neste trabalho. Microfluidica
envolve não só o fluido, mas também o trabalho de lidar com pequenos volumes, pequenos
tamanhos, baixo consumo de energia ou propriedades do fluido distintas das que,
normalmente, vemos nos fluidos macroscópicos.
A miniaturização leva a vários benefícios. Além dos benefícios encontrados quando se
miniaturiza componentes eletrônicos, como redução da perda de produtos, do custo e de
energia consumida, existem também: facilidade do controle de temperatura, uma vez que
pequenos volumes facilitam a difusão de calor; rapidez e precisão em reações químicas
devido à grande superfície de contato se comparado ao pequeno volume; e a utilização do
fluxo laminar contínuo em microcanais viabilizando novos métodos para realização de
procedimentos químicos [14].
Estes novos comportamentos em micro escala se devem ao baixo número de Reynolds
(Equação 1) que tende ao predomínio dos efeitos viscosos na dinâmica do fluido. Esta
propriedade é a mais importante para a determinação do regime de escoamento do fluido.
Nos casos em que o número de Reynolds é grande, a camada limite se desprende mais
facilmente fazendo com que o fluxo turbulento predomine, gerando flutuações caóticas nas
velocidades das partículas e aumentando a queda de pressão se comparado ao fluxo
laminar.
Equação 1 - Equação do número de Reynolds leva em consideração a densidade , a velocidade V, o comprimento
característico do canal L e da viscosidade do fluido .
A seguir, exemplos de fluxos predominantemente laminares em comparação a um
fluxo turbulento.
7
Figura 1 - Numero de Raynolds variando de 1 a 100 gerando fluxos laminares.
4
Figura 2 - Numero de Re = 10 .
Para um fluxo em um tubo de seção circular o número de Reynolds crítico, limiar ao
regime laminar para o turbulento, é de aproximadamente 2.300 [15].
Figura 3 – Simulação de escoamento laminar.
Figura 4 – Simulação de escoamento turbulento.
Fonte: NVIDIA. [7]
Considerando estas propriedades, a previsão do comportamento de fluidos em regime
laminar é muito mais segura do que em regime turbulento. Logo, simulações de dispositivos
8
microfluidicos são ferramentas confiáveis para desenvolvimentos e otimizações de novos
instrumentos tecnológicos.
No entanto, o trabalho com pequenas dimensões não tem somente vantagens. Esse
comportamento bem determinado pode ser considerado somente para sistemas
monofásicos, e suas simulações perdem confiabilidade à medida que o sistema passa a ser
mais complexo. A união de diferentes fases leva a dinâmicas ainda mais complexas como os
vórtices oriundos de forças centrífugas, chamados vórtices de Dean. Estes vórtices são muito
comuns em canais complexos, como microreatores químicos (Figura 5), mas também
aparecem em sistemas simples, como escoamentos de gotas por dentro de microcanais
(Figura 6).
Figura 5 – Exemplo de um microreator.
Figura 6 – Uma gota de água passando por dentro de um
tubo
Fonte: Chemical Engeneering Science. [16]
Fonte: Dynamics of microfluidic droplets. [17]
Portanto, a pesar de sistemas microfluidicos serem bem mais comportados que
sistemas turbulentos, os simuladores de microdispositivos têm um longo trabalho para
abranger todas essas dinâmicas. Um dos métodos de simulação que vem sendo muito
estudado e aperfeiçoado é o método DPD [18].
2.2. Método de Partículas Fluidicas (MPF)
Há diversos tipos de métodos de microsimulação [9], sejam eles baseados em malhas
computacionais, como o método CFD (Computational Fluid Dynamics) [19] [20], métodos
baseados em superfícies, como o método meshless Volume of Fluid [21], métodos híbridos,
como o SPH (Smoothed Particle Hydrodynamics) [22], ou mesmo métodos baseados em
partículas, como SDPD (Smoothed Dissipative Particle Dynamics) [23] e MD (Molecular
Dynamics) [24].
9
Figura 7 - Simulação de recirculação
em CFD.
Figura 8 - Simulação da
hidrodinâmica do peixe com
método meshless.
Fonte: Biomech Model Mechanobiol
[20]
Fonte: 16º POSMEC [21]
Figura 9 - Simulação SPH da quebra de
uma barragem.
Fonte: Tese Doutorado PUC/RJ [22]
Figura 10 – Simulação de um polímero em SDPD.
Fonte: IOPScience. [23]
Figura 11 - Membrana química simulada por MD.
Fonte: Journal of Physical Chemistry [24].
Cada um desses métodos de simulação possui diferentes nichos de abrangência.
Alguns possuem melhores resultados para simulações em larga escala, como moldagem de
metal líquido e dinâmica de ondas nas praias, outros modelam melhor iteração entre
partículas atômicas e propriedades físicas predominantes no ambiente microscópico.
De acordo com [18], há duas abordagens a se considerar para saber qual método
utilizar. A abordagem macro-microscópico e a abordagem micro-macroscópico.
Na primeira abordagem o estudo inicia-se com as equações de Navier-Stokes (Equação
2) para simulações de fluidos com grandes escalas de viscosidade, de densidade e de tempo.
A medida com que essas propriedades diminuem, outros métodos devem ser considerados,
passando por CFD, SPH e SDPD.
Equação 2 - Equação classica de Navier-Stokes para fluidos incompressíveis, na qual ρ é a densidade do fluido e η a
viscosidade, v o vetor velocidade, g o vetor aceleração de campo e P a pressão no fluido.
10
Na segunda abordagem as simulações tentam determinar grosseiramente a dinâmica
detalhada do sistema e conforme ele aumenta em tamanho, propriedades microscópicas
são retiradas por não serem tão influenciáveis no sistema, e assim os métodos vão se
alterando do MD às equações de Langevin (Equação 3) [25].
Equação 3 - Equação de Langevin, na qual, além da força externa Fi, o coeficiente de atrito e forças aleatórias R
influenciam na aceleração da partícula.
O método DPD proposto por (Steiner, Thomas (2009) [18]) leva em consideração essas
duas abordagem e indica um método mais robusto da representação realística de problemas
microfluidicos.
2.2.1. Método DPD
O método DPD que foi originado por Hoogerbrugge e Koelman tinha muitas limitações
como a não representação de fluidos complexos como polímeros e fluidos coloidais. O
método original foi baseado no método LGA (Lattice Gas Automata) [26] e em poucos anos
teve várias aplicações em mecânica estatística.
2.2.1.1. Equações de Movimento
A dinâmica do fluido é representada com iterações entre pares de partículas fluidicas.
Esta dinâmica é regida pela 2ª Lei de Newton.
̈
Equação 4 - Equação de movimento de Newton.
A força associada ao movimento das partículas pode ser subdividida em quatro forças
distintas.
∑
Equação 5 – Forças presentes na simulação do método DPD.
11
A primeira é a força conservativa (Equação 6), atribuída ao potencial suave que pode
ser aproximado pelo efeito da pressão entre diferentes partículas. Por ter esse potencial
suave, pode se considerar um tempo de integração grande no momento da simulação,
muito maior que no caso de MD.
{
(
)
Equação 6 - Força conservativa da partícula fluidica.
Esta força atua até certa distância ente as partículas, chamada de raio de corte (rc). O
parâmetro
define a magnitude da força entre partículas de diferentes grupos. A função
determina a influência do raio de corte entre as partículas e o vetor
a
direção normal da colisão entre as partículas. A figura a seguir explicita a atuação da força
conservativa em magnitude e direção.
Figura 12 - a) Função de magnitude da força conservativa dependente da distância entre as partículas. b) Exemplificação
de uma colisão entre partículas fluidicas e suas respectivas velocidades no momento do choque.
Fonte: Microfluid Nanofluid (2009) [27].
A segunda é a força dissipativa (Equação 7), oriunda das forças de atrito entre as
partículas, tem como objetivo diminui a quantidade de movimento relativo das partículas.
12
(
{
)
Equação 7 - Força dissipativa.
Nesta equação tem-se que
é o coeficiente de amortecimento que governa o
resfriamento do sistema. O produto escalar
é a componente da velocidade
relativa da partícula na direção da colisão, melhor visualizada na Figura 12. E a função de
ponderação
será descrita mais adiante, pois depende do teorema da flutuação-
dissipação que deve ser respeitado através das forças dissipativas e aleatórias.
A terceira é a força aleatória, que é estocástica e obedece à distribuição Normal [27]
tendendo a caotizar o sistema. Esta força é responsável pelo movimento Browniano das
partículas, ou seja, o movimento aleatório determinante para o aquecimento do sistema.
{
(
)
Equação 8 - Força aleatória.
O coeficiente
determina a amplitude da variável aleatória
. Essas duas últimas
forças, FD e FR, influenciam diretamente na viscosidade do fluido e agem de forma a deixar o
sistema em uma temperatura T.
A última força é a força externa que pode ser exemplificada pela força da gravidade.
A variável aleatória de
é quem insere a distribuição Gaussiana na equação da força,
e desta forma faz com que a esperança da força aleatória entre duas partículas seja sempre
zero. Juntamente com a escolha de um valor constante para o coeficiente de amortecimento
, estas propriedades satisfazem o teorema de flutuação-dissipação que determina os
outros parâmetros das forças,
e
(
)
(
)
(
), na qual
é
a constante de Boltzmann.
Com estas assertivas pode-se dizer, por saber-se que o sistema obedece à 3ª Lei de
Newton e que o número e o momento das partículas permanecerão constantes, que o
13
sistema é hidrodinâmico [18]. Deste modo podem ser determinados os parâmetros físicos
como pressão e densidade.
2.2.1.2. Abrangendo rotações
Com o intuito de ampliar o âmbito de sistemas que podem ser simulados pelo método
DPD, este foi melhorado para que abrangesse também sistemas que dependem do
momento rotacional das partículas. Logo
passou a representar não somente uma força
no sentido da normal da colisão, mas também uma combinação desta com a força de
cisalhamento perpendicular à normal da colisão.
Figura 13 - a) Forças dissipativas normais. b) Forças dissipativas perpendiculares.
Da mesma forma como na equação do movimento linear, o movimento angular é
regido pela Lei de Newton-Euler.
∑
∑
Equação 9 - Equação de euler do movimento de rotação, sendo que
partículas de diferentes tamanhos e
é um fator referente ao contato de
é a inércia da partícula fluidica.
Com o incremento desta componente na força dissipativa, a força aleatória deve se
modificar para que o teorema de flutuação-dissipação ainda valha e para que o sistema se
mantenha na mesma temperatura. Portanto, as novas equações das forças dissipativas e
aleatória são:
14
(
)(
)
(
)(
)
(
)
Equação 10 - Equação da força dissipativa com dinâmica rotacional.
̃
√
(
√
[
]
)
Equação 11 - Equação da força aleatória para suprir a nova força dissipativa.
Sendo que
e
são as componentes simétricas e assimétricas,
respectivamente, da matriz de ruído da distribuição Gaussiana.
O método DPD também foi aprimorado para tratar fluidos multifásicos e corpos
rígidos. O leitor encontrará maiores detalhes nas referências bibliográficas.
2.3. Programação paralela e GPU Computing
Na ultima década, as empresas fabricantes de processadores mudaram o conceito no
desenvolvimento de seus produtos. Pode-se observar que as gerações recentes de
processadores não possuem maiores frequência de operação que as gerações passadas, mas
sim, maior número de cores independentes. Esta mudança de pensamento foi oriunda de
dois motivos: o exagerado aquecimento e consumo de energia necessários para a
implementação de um novo processador com frequência ainda maior; e o surgimento
promissor dos processadores gráfico de propósito geral, as GPGPUs, com dezenas ou
centenas de processadores independentes.
Os processadores gráficos primários eram destinados à simples transferência de dados
para o monitor, avançando, nas versões posteriores, para o processamento de
transformações geométricas e iluminação. A revolução se iniciou na terceira geração de
GPUs, que com o processamento paralelo de texturas proporcionou um método rudimentar
de programação em suas arquiteturas paralelas [28]. Nas arquiteturas mais novas, as GPUs
evoluíram para unidades de processamento gráfico de propósito geral, as quais ganharam
linguagens de programação de alto nível [29] e algumas famílias de GPUs foram
desenvolvidas exclusivamente para este motivo, como é o caso da família TESLA [30].
15
Figura 14 - Evolução desde 2003 da largura de banda das
memórias comparando-se CPUs e GPUs
Figura 15 - Evolução desde 2003 da taxa de processamento
de diferentes tipos de GPU comparada à de CPUs.
Fonte: NVIDIA. [30]
Fonte: NVIDIA. [30]
O surgimento das linguagens de programação, que ainda estão em desenvolvimento,
como é o caso do CUDA e do OpenCL, ocasionou um grande impacto na comunidade de
desenvolvimento de softwares [31]. Criando-se o conceito de GPU Computing, que é o
emprego de placas gráficas modernas para processamento de aplicações altamente
paralelizáveis e com grande volume de dados [28].
Figura 16 – Separação das linguagens para GPUs em diferentes níveis de programação.
Fonte: GPU Programming Strategies. [29]
16
Além das linguagens, várias técnicas de programação estão em estudo para tornar
estas ferramentas ainda mais potentes. O OpenCL é uma linguagem open source que foi
criada recentemente pelo Khronos Group [32], para ser um ambiente de programação de
processamento paralelo voltado à sistemas híbridos, sejam eles constituídos de CPUs, GPUs
ou ambos. A linguagem CUDA foi criada pela NVIDIA [33], e é a mais estabelecida no
mercado por possuir vários recursos e alta performance, mas é exclusiva de suas próprias
GPUs.
A programação paralela está expandindo sua abrangência para todas as áreas da
computação, desde o projeto de hardwares, como os novos chips da INTEL, até em sites com
aceleração via GPGPU devido às linguagens como WebGL e WebCL. Logo, o paralelismo de
informações ou de tarefas tornou-se o futuro da computação [34], como é o caso do
desenvolvimento do Hoomd-blue [35] e outras aplicações em simulação de partículas [36].
2.4. Medição de Propriedades em Simulações Computacionais
A utilização de simulações computacionais normalmente está relacionada com
soluções numéricas de problemas dificilmente resolvidos analiticamente. Para que a
resposta vinda do computador seja estável e de acordo com a realidade é preciso que os
erros provenientes da simulação sejam minorados.
Esses erros são oriundos da discretização, no tempo e/ou no espaço, do problema e do
erro numérico devido à limitação na representação dos números inerente a todo
computador.
Visando contornar esse problema, há duas abordagens que podem ser executadas.
Quando se consegue um modelo matemático relativamente fiel, o número de subdivisões
necessárias no domínio pode ser reduzido deixando a cargo do modelo a correta
representação do sistema, ou então, pode-se subdividir o sistema tentando abranger até a
menor física relevante e deixar o modelo como sendo mais simplista, fazendo com que o
poder computacional seja o limitante na evolução do sistema. A abordagem dos programas
de simulação tem partido para a segunda opção levando a um nível microscópico o cálculo
da física das partículas.
17
No entanto, as propriedades físicas do fluido que se deseja obter, como a temperatura
e densidade, não são definidas na escala da simulação, mas sim em uma escala
macroscópica. Para avaliar essas propriedades são utilizados métodos estatísticos que a
partir dos estados de cada partícula determinam o estado macroscópico do sistema. Esse
problema é abrangido pela Teoria Cinética dos Gases.
A Teoria Cinética dos Gases iniciou-se em 1738 com a definição de que a pressão de
um gás ideal se deve as colisões dadas pelas partículas que o constituem na superfície do
volume que o abriga.
Assim, entre muitas propriedades, como energia cinética, temperatura, momento
vetor e outros, a distribuição da velocidade das partículas em um ambiente estático pode ser
determinada pela Lei de Distribuição das Velocidades Moleculares (distribuição de
velocidade de Maxwell-Boltzman). Dada pela Equação 12.
√
[
]
Equação 12 – Distribuição de velocidades em uma única direção.
Isso mostra que a média das velocidades devido à agitação das moléculas, em uma
determinada temperatura, é zero. Ou seja, a velocidade das partículas é puramente
aleatória e seu desvio padrão, nas três dimensões, pode ser determinado pela Equação 13.
√
Equação 13 – Desvio padrão de uma determinada velocidade nas partículas submetidas a uma temperatura.
2.4.1. Filtro de Velocidades
Para um sistema cujas partículas possuem velocidade relativamente baixa, tem-se a
necessidade de reduzir esse fenômeno de aleatoriedade da energia cinética das partículas
devido à temperatura. Neste caso é possível utilizar um filtro, igualmente feito em filtragem
18
de imagens, o qual a velocidade do fluido em certo ponto pode ser obtida pela média das
velocidades das partículas ao redor.
Esta técnica é utilizada em filtragem de imagens para remover ruídos pontuais da
imagem. Como mostrado na figura a seguir.
Figura 17 – Imagem com ruidos chamados “Salt and
Pepper” Noise.
Figura 18 – Imagem após a aplicação do filtro de média.
Fonte: Lode's Computer Graphics Tutorial - Image Filtering.
[37]
Fonte: Lode's Computer Graphics Tutorial - Image Filtering.
[37]
Este filtro é realizado pela convolução dos dados de entrada com uma função, que no
caso é a função de média. Esta operação é sistematizada da seguinte forma:
Equação 14 – Operação de convolução, na qual h é a função de saída, f a função de entrada e g a função filtro.
Figura 19 – Sistema de blocos da convolução.
Fonte: Machine Vision – Capitulo 4, pg 116. [38]
A convolução de duas funções contínuas é determinada pela seguinte fórmula.
19
∫ ∫
Equação 15 – Função de convolução de funções contínuas.
No entanto, em imagens é mais comum utilizar a formulação discreta, já que o espaço
é discretizado em pixels.
[
]
[
]
[
]
∑∑ [
] [
]
Equação 16 – Equação discreta da convolução.
Exemplificando, admita-se uma imagem com ruídos. Cada pixel (P) da imagem será
convoluido com uma matriz de convolução (A) para a eliminação do ruído. Esta operação é
realizada da seguinte forma.
Figura 20 – Exemplo de convolução discreta em figuras.
Fonte: Machine Vision – Capitulo 4, pg 117. [38]
A imagem mostra o pixel (h[i,j]) sendo atualizado pela convolução de sua matriz P
(f[i,j]) pela matriz de convolução A (g[i,j]).
20
Neste caso temos uma acumulo de brilho (ou energia cinética no caso do estudo das
velocidades, i.e. no caso aqui abordado) na utilização deste método. Portanto, ao atualizar
uma posição h[i,j] é necessário normalizar o resultado. Esta normalização é realizada
dividindo-se o resultado (h[i,j]) pela soma das grandezas de g[i,j].
[
]
[
(a)
]
(b)
Figura 21 – Exemplo de matriz de convolução normalizadas. a) uma matriz de convolução de 3x3 e em b) uma matriz de
convolução 5x5.
É importante ressaltar que as matrizes de convolução devem ter uma célula central, a
qual influenciará diretamente no pixel examinado da imagem. Ou seja, as matrizes precisam
ser de dimensões ímpares.
2.5. Hoomd-blue
Hoomd-blue é um software open source, desenvolvido em CUDA C que usa
processadores gráficos (GPU) da NVIDIA para simulação de dinâmica de partículas.
Desenvolvido pelo departamento de engenharia química da universidade de Michigan e teve
sua primeira versão em 2009.
Desde lá, o programa obteve muitas contribuições e com isso cresceu rapidamente.
Atualmente, na versão 0.10.1 (momento em que este trabalho foi iniciado), possui
implementado, diversos métodos de partículas, várias condições de contorno, vários
métodos de integração, além de ser facilmente configurável por arquivos externos ao código
e exportar resultados tanto para análise quanto para a visualização da simulação. A seguir,
imagens dos sistemas que podem ser simulados no software.
21
Figura 22 - Slide de demonstração das possíveis dinâmicas que podem ser simuladas pelo Hoomd-blue.
Fonte: Slides Advancing GPU Molecular Dynamics. [39]
Alguns dos métodos de partículas contidos no software são: CGCMM, EAM (embedded
atom mthod), Gaussian, Lennard-Jones, Yukawa, além do DPD que será estudado neste
trabalho.
O projeto Hoomd-blue, tem a finalidade de facilitar o contato com a tecnologia da
simulação de dinâmica de partículas utilizando-se de processamento paralelo. Deste modo, a
comunidade interessada está alavancando suas pesquisas e desenvolvimentos, devido a sua
versatilidade na gama de opções de métodos para simulação, a sua velocidade com o uso do
processamento paralelo em GPUs e CPUs, e a sua facilidade, com modos de configurações
práticos.
Apesar das facilidades na execução de um sistema desejado, devido ao seu método
básico de configuração via Python script, o software também possui recursos avançados.
Com isso, tem-se a possibilidade de: implementação de plugins para uso de programação
adicional; utilização do software VMD para visualização da simulação online ou gravada em
arquivo XML; exportação de dados das posições das partículas para arquivo DCD para
posterior análise; além de outros recursos.
22
2.5.1. Método DPD no Hoomd-blue
A programação de um método de partículas é, basicamente, dividida em três etapas:
1) dada uma partícula do sistema, obter a lista de suas partículas vizinhas; 2) calcular as
forças oriundas da interação entre essas partículas; 3) realizar a dinâmica da força aplicada
às partículas através do tempo discretizado por meio de um método de integração.
A evolução do método DPD, no Hoomd-blue, foi concretizada com base no estudo
realizado por Phillips et al [40]. Neste estudo, o autor indica uma forma mais rápida,
utilizando GPUs, de obter números estocásticos para o cálculo das forças aleatórias.
Há três maneiras de simular o comportamento de partículas fluidicas, no software,
pelo método DPD. Uma forma é o método normal, semelhante ao mostrado na seção 2.2.1,
outro é excluindo-se a assertiva da natureza térmica e o ultimo, substituindo-se o cálculo da
força conservativa por um calculo que inclui a dinâmica entre partículas de Lennard Jones.
2.5.1.1. DPD Termostato (clássico)
Análogo ao método descrito por Steiner, o programa implementado utiliza-se das
equações de força conservativa, força dissipativa e força aleatória para o cálculo das forças
de contato entre duas partículas. A seguir, as equações exatas utilizadas pelo simulador
Hoomd-blue no método DPD.
(
)
(
(
(
(
(
)
)
(
)
)( ̂
√ √
)
{
(
)
)
(
)
)
Equação 17 - Equação do método DPD para cálculo da força entre partículas do software Hoomd-blue.
Pode-se notar a semelhança entre as equações. Considerando-se que o teorema de
flutuação-dissipação já está incluso, as diferenças estão no
respectivamente nomeados anteriormente de
e no
que foram
e .
23
No script de configuração, a função a ser utilizada é pair.dpd cujos parâmetros são o
raio de corte
, a temperatura do sistema , a magnitude da força conservativa
constante de amortecimento
, podendo
e
ea
serem diferentes para cada grupo de
partículas criado.
Esta função deve ser utilizada conjuntamente com o método de integração via
Velocity-Verlet, integrate.nve, que já abrange as necessidades termodinâmicas do método. A
utilização de outro método de integração leva o sistema à instabilidade.
2.5.1.2. DPD Conservativo
Este método se diferencia pela não utilização dos fatores oriundos da dinâmica
térmica para o cálculo da força entre as partículas. Portanto, os fatores
e
não são
utilizados.
A força conservativa é então calculada da seguinte maneira:
{
Equação 18 - Equação da força conservativa deste método DPD.
Deste modo, o sistema deixa as propriedades microscópicas de lado e passa a se
comportar com uma dinâmica macroscópica.
2.5.1.3. DPD com Lennard Jones
Neste método, a termodinâmica do sistema retorna, mas a força conservativa não é
dada pelo gráfico da Figura 12a, mas sim, dependente da função de Lennard Jones clássica,
para o potencial dessas interações entre partículas.
Tem-se então a seguinte alteração na equação de
{
[( )
:
( ) ]
Equação 19 - Equação da força conservativa utilizando potencial de Lennard Jones.
24
A implementação do método DPD utilizando-se do potencial de Lennard Jones, possui
uma diferente abordagem que facilita a parametrização do sistema de partículas que será
simulado [41].
Portanto, estas opções tornam o Hoomd-blue ainda mais abrangente, no entanto,
neste trabalho será utilizado o método DPD clássico.
2.6. Modelo matemático do perfil de velocidades
A hidrodinâmica proposta por Hagen-Poiseuille modela muito bem um fluxo
estacionário de um líquido no interior de um microcanal cilíndrico. Uma vez, que o modelo
depende de suposições garantidas por um microsistema, como, fluxo laminar, viscoso e
incompressível visto na seção 2.1.
A equação/lei proposta por eles é derivada da equação de Navier-Stokes, a qual em
um microtubo os termos não lineares são desprezados e o perfil de velocidades depende do
gradiente de pressão em uma lâmina frontal do líquido [42]. O perfil de velocidades é
oriundo da diferença de pressão entre os dois lados do tubo, e permanece parabólico,
devido às forças de atrito e à condição de não deslizamento nas paredes do tubo [43].
Figura 23 – Perfil de velocidade determinado pela equação de Hagen-Poiseuille para um microcanal.
Fonte: Journal of the American Chemical Society [43].
O perfil parabólico, para um fluido viscoso, incompressível e com escoamento
permanente, é determinado pela equação de Hagen-Poiseuille descrita a seguir:
Equação 20 - Equação de Hagen-Poiseuille do perfil de velocidade do fluxo em um microcanal. O qual
do líquido, o tamanho do tubo, e metade da espessura do tubo.
é a viscosidade
25
Esta equação foi utilizada por Gösch et. al. [44] para obter o perfil de velocidade em
um microcanal de seção retangular de 50x50 μm². Ele lembra que, deve-se garantir um
número de Reynolds (Equação 1) abaixo de 2000, e para isso foi admitida uma velocidade
máxima de V = 70 mm/s.
2.6.1. Diâmetro Hidráulico
Os experimentos são realizados comumente em tubos de seção cilíndrica, pois são os
mais utilizados, no entanto, no caso aqui descrito vamos utilizar um tubo de seção quadrada.
Para realizar a correlação entre um tubo de seção quadrada e um tubo cilíndrico
utiliza-se o diâmetro hidráulico para os cálculos das formulas que levam em consideração o
raio do tubo.
Deste modo tem-se a correlação dos seguintes tipos de seção do tubo com seus
respectivos diâmetros hidráulicos na Figura 24.
Figura 24 - Correlação do diâmetro hidráulico com diversos tipos de seção de tubo.
Com isso fechamos o ambiente físico e matemático que serão necessários na validação
do experimento realizado.
26
3. Métodos
3.1. Metodologia de validação do software
O software será validado com base na comprovação de um dos experimentos
realizados por Steiner [18] chamado The Gravitational Method, o qual se utiliza da equação
de Hagen-Poiseuille para comparar o perfil de velocidade de um determinado sistema
microfluidico.
Este experimento foi definido por 100.000 partículas no interior de um tubo de seção
quadrada que sofrem a ação de uma força externa constante
, a qual determina o
gradiente de pressão do tubo gerando um perfil de velocidade. A estrutura do tubo é
constituída por partículas fixas, que geram a propriedade de não deslizamento necessário
para o método. A condição de contorno do tubo, no sentido do fluxo, é definida pelas
condições de Lees-Edwards, que periodizam o comprimento, fazendo com que as partículas
que saiam de um lado do tubo reentrem do outro lado. Deste modo, a conservação de
matéria, indispensável pelo método DPD, é alcançada.
3.2. Configuração do sistema proposto
O sistema, determinado para a validação da metodologia, será interpretado pelo
software Hoomd-blue, através de seu arquivo básico de configuração em linguagem Python
(Anexo A – Arquivo de Configuração do Hoomd-blue).
O sistema de partículas será composto por 3 tipos de partículas: as partículas do tipo
A, serão o fluido com a dinâmica regida pelo método DPD; as partículas do tipo B, serão fixas
e determinarão as paredes do tubo; e as partículas do tipo C, serão uma fatia do fluido para
melhor visualização do efeito de comprimento de entrada e ter uma noção do perfil de
velocidades.
O tubo quadrado terá tamanho de 13x13 partículas e comprimento de 39 partículas,
sendo que as partículas possuirão raio de 1 unidade de comprimento. Todas as partículas
serão dispostas uniformemente em um paralelepípedo cujas faces, superior e lateral, serão
do tipo B, e as do interior como partículas do tipo A, exceto, a fatia de partículas mais
próxima da face periódica que serão do tipo C.
27
A condição de contorno periódica já é, por default, realizada pelo software e não terá a
necessidade de ser configurada. As dinâmicas das partículas A e C serão geridas pelo método
DPD clássico com
e
. O coeficiente de amortecimento
entre as partículas deve
ser determinado pela densidade e temperatura do sistema de acordo com a expressão de
flutuação-dissipação
, mas no caso será adotado
por não ser necessária
a correta representação de algum fluido específico e para que a simulação seja similar à
apresentada na literatura. O parâmetro de repulsão
e determinou um número relevante de
foi estudado com cautela por Steiner
que minimiza o erro relativo da viscosidade
do fluido simulado. A viscosidade do sistema pode ser obtida alterando-se os seguintes
parâmetros: amplitude da variável estocástica ; magnitude do parâmetro de repulsão entre
partículas
equação
; e a densidade
, que também determina a pressão do sistema através
. No caso do Hoomd-blue o parâmetro
será calculado da seguinte
forma:
∑
Equação 21 – Equação da densidade.
Sendo, no numerador, a somatória da massa de todas as partículas do tipo A e C, e, no
denominador, o volume do paralelepípedo em número de partículas incluindo-se a
porcentagem do espaçamento inicial entre elas.
O experimento de Steiner determina uma densidade de 5 para suas simulações. Deste
modo, o sistema foi alterado atribuindo-se um valor de massa às partículas de modo que o
fluido chegasse à mesma densidade.
O tempo será discretizado em passos de 0,002 unidades de tempo, que foi constatado
ser suficiente para se obter um preciso controle da temperatura do sistema. Inicialmente,
todas as partículas serão iniciadas automaticamente com velocidades estocásticas para
determinação da temperatura desejada.
28
Igualmente realizado no trabalho de validação, serão executadas quatro simulações.
Todas elas terão alteração somente na força de campo aplicada no fluido. Estas forças serão:
0,0253; 0,0355; 0,0480; 0,0600 [unidades de força].
Todas as simulações serão configuradas para executarem por 40.000 iterações.
Os dados das posições das partículas, a cada iteração, serão salvos pelo próprio
Hoomd-blue em arquivos “XML”, que serão trabalhados via MatLab para obtenção do perfil
de velocidade (Anexo B – Programa em MatLab).
3.3. Tratamento dos dados
Os dados obtidos de velocidade das partículas em uma fatia da seção do tubo serão
comparados com os resultados simulados e teóricos da validação da literatura.
No entanto, para melhor visualização, os dados serão filtrados utilizando-se a
convolução discretizada da velocidade das partículas com uma máscara de média para
minimizar o ruído gaussiano das velocidades devido à temperatura do sistema, como
indicado pela seção 2.4.
Como no método utilizado em imagens, os pontos utilizados para a descoberta da
velocidade do fluido serão os centros das células de uma malha alocada na seção transversal
do tubo a uma distância d do início.
Figura 25 - Malha no interior do tubo cujas celulas informarão os valores médios da velocidade do fluido na região em
que está alocada.
Mas, a técnica que será utilizada não levará em conta outra matriz para a função de
convolução, mas sim realizará uma média das partículas com raio menor a uma distância
específica (raio de corte) do centro da célula a ser atualizada. Como ilustrado a seguir:
29
Figura 26 - Centro da célula com raio de corte definido.
Com o raio de corte definido, a velocidade do fluido na região da célula em questão
será calculada pela média das velocidades das partículas no interior do espaço determinado
pelo raio de corte. Desta forma, as partículas em vermelho serão consideradas no cálculo e
as partículas em verde serão desconsideradas.
Após a filtragem dos dados, a matriz das velocidades médias sofrerá uma interpolação
dos dados para minorar os problemas provindos da discretização em partículas do sistema.
Assim, os gráficos obtidos do perfil de velocidades serão comparados com os gráficos
presentes no trabalho de Steiner.
3.3.1. Programa de Tratamento
O sistema de tratamento de dados foi desenvolvido em MatLab, por ela ser uma
ferramenta própria para trabalhar com grande número de dados encadeados em forma de
tabela e por ter uma grande variedade de gráficos que melhor representam o problema.
Assim, a lógica utilizada pode ser dividida em quatro etapas. São elas: Importação,
Redução, Convolução e Exibição.
A etapa de importação resume-se a importar os arquivos XML gerados pelo programa
Hoomd-blue de forma a facilitar a manipulação dos dados no interior do programa de
tratamento.
Esta etapa foi basicamente realizada por pesquisa bibliográfica, já que arquivos XML
são utilizados em várias aplicações como arquivo padrão de disposição de dados. Assim,
utilizou-se o programa xml2struct desenvolvido por Wouter Falkena em 2010 [45].
30
Com os dados corretamente dispostos, o programa de determinação do perfil de
velocidades recebe os seguintes parâmetros.
 Dados do XML já tratados.
 Distância da origem na qual será inserida a malha de verificação.
 Raio de corte.
 Porcentagem de refinamento da malha se comparado à área total da
seção transversal do tubo.
Visando diminuir o número de partículas a serem verificados, a etapa de Redução
realiza duas eliminações na tabela de dados. Uma é no início do programa, visando diminuir
o sistema de partículas para as partículas que estão a uma distância rC (raio de corte) da face
onde a malha foi alocada.
Figura 27 - Primeira redução.
A outra redução é realizada em cada elemento da malha. Que ao iniciar o
processamento do elemento, o sistema de partículas é novamente reduzido às partículas
que estão somente a uma distância rC do centro da célula
Figura 28 - Segunda redução.
31
Com cada célula e seus respectivos sistemas locais de partículas, a etapa de
Convolução é realizada. Esta etapa é similar à convolução discreta em processamento de
imagens com uma malha de filtro mediano. No entanto, a malha é de tamanho alterável
conforme o número de partículas na proximidade da célula em execução.
Caso haja uma partícula no interior da região, a velocidade na direção do fluxo é
somada na “velocidade” da célula. Quando não houver mais partículas ao redor daquela
célula, o montante da velocidade é dividido pelo número de partículas que foram
consideradas. Esta operação visa normalizar a quantidade de movimento da operação de
convolução, conforme descrito na seção 2.4.1.
A última etapa é a de Exibição dos dados. Foram escolhidos dois gráficos para serem
exibidos que retratam o perfil de velocidade calculado pela malha. Um gráfico demonstra o
perfil de velocidades em três dimensões averiguando a sua forma geral, e o outro gráfico de
duas dimensões exibe o perfil de velocidade em uma linha transversal da malha.
Figura 29 - Seção do tubo de partículas denotando a linha transversal utilizada para obtenção do perfil de velocidade 2D.
Este último gráfico é o que será comparado com os gráficos determinados por Steiner.
4. Resultados e Discussões
4.1. Resultado qualitativo
A execução da simulação foi realizada pelo software em 10 segundos em um
computador com GPU NVIDIA GTX 560 Ti. As posições das partículas a cada período de 100
32
iterações foram exportadas também para um arquivo “DCD” para visualização pelo software
VMD [46].
A seguir imagens da execução do resultado do sistema pelo software VMD.
Figura 30 – Ilustração do tubo retangular em roza,
partículas do tipo B, e o fluido em azul, partículas do tipo
A.
Figura 31 – Perfil do tubo com partículas fluidicas com
fluxo na direção Z.
Figura 32 – O sistema com as paredes translucidas em sua
posição inicial. As partículas em roxo são do tipo C.
Figura 33 – Sistema após 854 frames mostrando pela
posição das partículas do tipo C o fluxo parabólico.
33
Figura 34 – No mesmo momento que a figura anterior, agora em perspectiva 3D, mostrando o fluxo parabólico em todas
as faces do tubo.
As duas primeiras imagens exibem a configuração do tubo e das partículas em seu
interior. A terceira indica a posição inicial das partículas para comprovar sua evolução com o
tempo. As últimas mostram o perfil da posição das partículas no momento em que melhor
despontam o perfil de velocidade, denotando o perfil da região de entrada do tubo.
Este perfil pode ser mais bem visualizado na seguinte sequência de figuras. As
seguintes simulações já foram realizadas com a configuração indicada por Steiner. Os dados
coletados dessas simulações serão as utilizadas na validação do software.
Figura 35 - Sequência de formação do perfil de entrada na simulação com força de campo de 0,0253.
34
Figura 36 - Sequência de formação do perfil de entrada na simulação com força de campo de 0,0355.
Figura 37 - Sequência de formação do perfil de entrada na simulação com força de campo de 0,0480.
Figura 38 - Sequência de formação do perfil de entrada na simulação com força de campo de 0,0600.
4.2. Resultado quantitativo
Realizados os experimentos finais na configuração de Steiner, os dados foram
angariados e tratados pelo software MatLab.
A seguir os dados da literatura utilizando 4 campos de força diferentes.
35
Figura 39 - Perfil de velocidades obtidos por Steiner, pg 46 [18].
Este gráfico retrata uma das experiências realizadas por Steiner, que é a validação pelo
método gravitacional. Pode-se verificar que foram utilizadas quatro diferentes forças de
campo, sendo que, com o aumento das forças, há uma maior variação no perfil de
velocidades.
Os pontos triangulares do gráfico são as amostras de sua experiência, e a linha
contínua é a resposta analítica obtida pela função de Hagen-Poiseuille.
A seguir, os quatro gráficos deste experimento gerados pelo MatLab, um para cada
força de campo, exibindo os dados tratados provenientes das simulações realizadas no
Hoomd-blue.
36
Figura 40 - Resultado obtido com força de campo de 0,0253 unid. de força.
Estes são os dois gráficos obtidos na simulação com força de campo de 0,0253
unidades de força.
Inicialmente no gráfico em três dimensões é possível verificar o paraboloide formado
pela malha de velocidades indicando uma correta representação do perfil de velocidades. Na
parte abaixo do gráfico há curvas de nível que indicam um ponto de máximo no centro da
figura.
Verifica-se também que ao aproximar das bordas, as curvas de nível passam de
círculos para retângulos de cantos arredondados. Isso denota o diâmetro hidráulico do tubo
de seção retangular.
Um erro que pode ser averiguado é o pico de velocidade no ponto superior direito da
figura. Este pico pode ser proveniente da interpolação dos dados próximo à borda, gerando
outro ponto de inflexão. Mas também pode ser da própria simulação, que por ter poucas
partículas na seção transversal do tubo e pela baixa velocidade constatada, um ruído pode
ter interferido de forma significativa na interpolação.
Já o gráfico de duas dimensões tem perfil parabólico e velocidade máxima de 0,27
unidades de velocidade, ficando muito parecida com a velocidade máxima obtida, tanto no
experimento de Steiner, quanto na resposta analítica.
37
Pode-se notar também que nos cantos do perfil, as velocidades das partículas se
tornam negativas. Ou seja, há uma recirculação do fluido no local próximo à borda.
Este fenômeno é possível e é encontrado frequentemente em estudos de camada
limite, como pode ser observado na seguinte figura:
Figura 41 - Fluxo reverso, indicando velocidade negativa na proximidade da borda.
O próximo gráfico expõe os dados da simulação com força de campo de 0,0355
unidades de força.
Figura 42 - Resultado obtido com força de campo de 0,0355 unid. de força.
O gráfico em três dimensões tem as mesmas propriedades já comentadas da
simulação anterior, no entanto, não há mais o erro de velocidade na borda do tubo. Isso
pode ter ocorrido pois neste ponto a velocidade já é mais aparente.
38
Também é possível verificar uma acentuação nas curvas de nível que ficaram ainda
mais retangulares próximas à parede.
O gráfico em duas dimensões do perfil de velocidade também tem o perfil parabólico e
seu pico tem uma velocidade de 0,42. Neste caso a velocidade ficou um pouco distante do
obtido por Steiner, no entanto, o aumento da velocidade máxima com o aumento da força
de campo foi respeitada.
Observa-se, também, uma lateralização do valor máximo, o que ocasionou a diferença
das velocidades nas duas bordas. Uma possibilidade deste acontecimento, é um acúmulo de
partículas do lado esquerdo melhorando a estimativa da velocidade deste lado, em
detrimento à estimativa do outro lado.
A seguir o gráfico da simulação com força de campo de 0,048 unidades de força.
Figura 43 - Resultado obtido com força de campo de 0,0480 unid. de força.
Neste caso vemos uma simetria maior dos gráficos devido à alta velocidade relativa
das partículas do centro em relação às partículas periféricas.
No gráfico em três dimensões observa-se os mesmos pontos levantados no gráfico
anterior.
39
No gráfico em duas dimensões verifica-se a simetria tanto na centralização do ponto
de máximo da parábola, como das velocidades próximas à parede.
A velocidade máxima obtida foi de 0,51 unidades de velocidade retornando a
proximidade com os dados obtidos da literatura.
A seguir o gráfico da simulação realizada com força de campo de 0,06 unidades de
força.
Figura 44 - Resultado obtido com força de campo de 0,0600 unid. de força.
Pode-se verificar uma grande simetria nos dois gráficos. No entanto, ambos gráficos
possuem uma variação no lado esquerdo. Esta variação, nos dois gráficos, pode ter sido
resultada da mesma causa que também foi encontrada na primeira simulação. Esta causa é
muito provavelmente determinada pela escassez de partículas próximas aos cantos. Com
isso gera-se um erro maior na medição da velocidade que ocasiona esta inflexão presente
nestes dois gráficos.
O gráfico em duas dimensões, a pesar de sua pequena descentralização, confirma as
propriedades obtidas na simulação anterior.
A velocidade máxima obtida neste caso foi de 0,65 unidades de velocidade, ficando
razoavelmente próxima à velocidade obtida pela literatura.
40
É interessante também averiguar a seguinte figura da última simulação.
Figura 45 - Vista superior do grafico 3D da última simulação.
Pode-se evidenciar neste gráfico o diâmetro hidráulico do tubo e a transformação de
curvas de nível em círculos para retângulos.
5. Trabalhos Futuros
Um possível trabalho futuro é desenvolver um software próprio de simulação DPD
para ganhar know-how e contribuir para um software ainda maior que abranja todas as
simulações de partículas, deixando à caráter do usuário qual método escolher.
Além disso, o aprofundamento em técnicas de paralelização sejam elas por GPUs,
CPUs, hibridas ou em rede para tornar viável estes tipos de simulação alcançando patamares
em tempo real, independente de qual for o número de partículas no domínio.
6. Conclusão
Os resultados qualitativos atingiram os objetivos na execução do software Hoomd-blue
simulando um fluido escoando por dentro de um tubo de seção quadrada. A dinâmica do
fluido foi regida pelo método DPD, que após ajustes em seus parâmetros, salientou o perfil
de velocidade das partículas exibindo o perfil parabólico da posição das partículas roxas,
como o esperado.
41
A quantificação do perfil de velocidades foi satisfatória, sendo que o perfil parabólico
de velocidades foi respeitado em todas as simulações e a velocidade máxima obtida ficou
próxima à esperada observando os dados da literatura. Certamente o resultado das
simulações não culminaria em dados exatamente iguais ao da literatura, pois, o modo de
cálculo dos dois programas, considerando hardware e software, é diferente, além de que, a
aleatoriedade da velocidade determinada pela temperatura do sistema aumenta a
dificuldade de reprodução de igualdade de dados de duas simulações com os mesmos
parâmetros.
Em geral o estudo foi de grande proveito e os resultados obtidos foram os esperados.
Portanto, por meio deste trabalho, indico a utilização do software Hoomd-blue para a
validação de experimentos com simulação microfluidica.
42
7. Bibliografia
1.
RIBAS, R. P. Microssistemas Integrados (MEMS). [S.l.]: [s.n.].
2.
OJEA, R. B. Numerical and Experimental Analyses of Single And Two-Phase
Microfluidic Flows with Implications in Microreactors. Tarragona: [s.n.], 2011.
3.
CHE, Z.; NGUYEN, N.-T.; WONG, T. N. Analysis of chaotic mixing in plugs moving in
meandering microchannels. PHYSICAL REVIEW E 84, 066309, 2011.
4.
CHE, Z. et al. Numerical investigation of upstream pressure fluctuation during
growth and breakup of pendant drops. Chemical Engineering Science, 2011.
5.
CHOI, S.-W.; ZHANG, Y.; XIA, Y. A Temperature-Sensitive Drug Release System
Based on Phase-Change Materials. Angew. Chem. Int., 49, 7904 –7908, 2010.
6.
JÖNSSON, H. Microfluidics for lab-on-a-chip applications. [S.l.]: [s.n.].
7.
PFAFF, T. et al. Scalable Fluid Simulation using Anisotropic Turbulence Particles.
NVIDIA. [S.l.]. 2010.
8.
CLAVET, S.; BEAUDOIN, P.; POULIN, P. Particle-based Viscoelastic Fluid Simulation.
Eurographics/ACM SIGGRAPH Symposium on Computer Animation, 2005.
9.
MEAKIN, P.; XU, Z. Dissipative Particle Dynamics and Other Particle Methods for
Multiphase Fluid Flow in Fractured and Porous Media. 6th International Conference on
CFD in Oil & Gas, Metallurgical and Process Industries. [S.l.]: [s.n.]. 2008.
10.
OWENS, J. D. et al. A Survey of General-Purpose Computation on Graphics
Hardware. Computer Graphics Forum, 26(1):80–113, Março 2007.
11.
FUJIOKA, H.; TAKAYAMA, S.; GROTBERG, J. B. Unsteady propagation of a liquid plug
in a liquid-lined straight tube. PHYSICS OF FLUIDS 20, 062104, 2008.
12.
KREUTZERA, M. T. et al. Multiphase monolith reactors: Chemical reaction
engineering of segmented flowin microchannels. Chemical Engineering Science 60 5895
43
– 5916, 2005.
13.
FOX, R. W.; MCDONALD, A. T.; PRITCHARD, P. J. Introdução à Mecânica dos
Fluidos. 6ª. ed. [S.l.]: LTC, 2006.
14.
VAN DAM, R. M. Introduction to Microfluidics. In: VAN DAM, R. M. Solventresistant elastomeric microfluidic devices and applications. [S.l.]: California Institute of
Technology, 2006. Cap. 2.
15.
KARNIADAKIS, G. M.; BESKOK, A.; ALURU, N. Essentials of Fluidic Transport
Phenomena at Small Scales. In: KARNIADAKIS, G. M.; BESKOK, A.; ALURU, N. Microflows
and Nanoflows. [S.l.]: [s.n.], 2005. Cap. 2.
16.
BOTHE,
D.;
LOJEWSKI,
A.;
WARNECKE,
H.-J.
Fully
resolvednumericalsimulationofreactivemixinginaT-shaped micromixer using parabolized
species equations. Chemical Engineering Science.
17.
BAROUD, C. N.; GALLAIREB, F.; DANGLA, R. Dynamics of microfluidic droplets DOI:
10.1039/c001191f. Lab on a Chip, Janeiro 2010.
18.
STEINER, T. Dissipative Particle Dynamics - Simulation of Microfluidic Systems
With Fluid Particle Methods on High Performance Computers. [S.l.]: Universität Freiburg
- Departament of Microsystems Engineering, v. V, 2009.
19.
DEL FRARI, R. V. Malhas Computacionais para Simulação Numérica de
Escoamentos de Fluidos Entre Cilindros com Excentricidade. X Encontro Gaúcho de
Educação Matemática. Ijuí/RS: [s.n.]. 2009. p. 8.
20.
FENG, R. et al. Viscous flow simulation in a stenosis model using discrete particle
dynamics: a comparison between DPD and CFD. Biomech Model Mechanobiol 11:119–
129, 2012.
21.
ESTACIO, K. C.; FERREIRA, V. G.; NONATO, L. G. Um Método Meshless Para
Simulação de Escoamento de Fluidos em Cavidades de Moldes; Mangiavacchi,
44
Norberto. Uberlândia: Universidade Federal de Uberlândia, 2006.
22.
NETO, A. P. Uma abordagem lagrangeana para simulação de escoamentos de
fluidos viscoplásticos e multifásicos. Rio de Janeiro: [s.n.], 2007.
23.
LITVINOV, S.; HU, X. Y.; ADAMS, N. A. Numerical simulation of tethered DNA in
shear flow. JOURNAL OF PHYSICS: CONDENSED MATTER, 2011.
24.
LIN, J.; NOVAK, B.; MOLDOVAN, D. Molecular Dynamics Simulation Study of the
Effect of DMSO on Structural and Permeation Properties of DMPC Lipid Bilayers. The
Journal of Physical Chemistry, 2011.
25.
CARVALHO, J. D. C.; CANEPPELE, N.; FIGUEIREDO, C. M. Um novo método para
resolver a equação de Langevin aplicada à dispersão de poluentes atmosféricos em
regime de turbulência Gaussiana. ACTA Scientia, v. 5, n. 1, 2003.
26.
ZHENGA, Z. et al. A dual-scale lattice gas automata model for gas–solid two-phase
flow in bubbling fluidized beds. Computers and Mathematics with Applications, 2011.
27.
STEINER, T. et al. Simulation of advanced microfluidic systems with dissipative
particle dynamics. Microfluid Nanofluid 7:307–323, Janeiro 2009.
28.
IKEDA, P. A. Um estudo do uso eficiente de programas em placas gráficas. São
Paulo/BR: [s.n.], 2011.
29.
BRODTKORB, A. R.; HAGEN, T. R.; SÆTRA, M. L. GPU Programming Strategies and
Trends in GPU Computing. Journal of Parallel and Dritributed Programming, Oslo, 2012.
30.
NVIDIA. NVIDIA CUDA C - Programming Guide. [S.l.]: [s.n.], v. Ver. 4.2, 2012.
31.
KIRK, D.; HWU, W.-M. Programming Massively Parallel Processors: A Hands-on
Approach. 57p. Massachusetts: [s.n.], 2010.
32.
KHRONOS GROUP. Khronos Group - Connecting Software to Silicon. The Khronos
Group Inc., 2012. Disponivel em: <http://www.khronos.org>. Acesso em: jun. 2012.
45
33.
NVIDIA CORPORATION. Plataforma de Computação Paralela | CUDA | NVIDIA. Site
da
NVIDIA
Corporation,
2012.
Disponivel
em:
<http://www.nvidia.com.br/object/cuda_home_new_br.html>. Acesso em: jun. 2012.
34.
OWENS, J. D. et al. GPU Computing. [S.l.]: IEEE, 2008.
35.
THE GLOTZER GROUP. Highly Optimized Object-oriented Many-particle Dynamics
(Hoomd-blue). Hoomd-blue, 2011. Disponivel em: <http://codeblue.umich.edu/hoomdblue>. Acesso em: jun. 2012.
36.
ANDERSON, J. A.; LORENZ, C. D.; TRAVESSET, A. General purpose molecular
dynamics simulations fully implemented on graphics processing units. Journal of
Computational Physics 227(10): 5342-5359, Maio 2008.
37.
VANDEVENNE, L. Image Filtering. Lode's Computer Graphics Tutorial, 2004.
Disponivel em: <http://lodev.org/cgtutor/filtering.html>. Acesso em: 01 jun. 2013.
38.
JAIN, R.; KASTURI, R.; SCHUNCK, B. G. Image Filtering. In: JAIN, R.; KASTURI, R.;
SCHUNCK, B. G. Machine Vision. [S.l.]: McGraw-Hill, Inc., 1995. Cap. 4, p. 112-139.
39.
ANDERSON, J. A. et al. Advancing GPU Molecular Dynamics: Rigid Bodies in
Hoomd-blue. The Glotzer Group University of Michigan. [S.l.]. 2011.
40.
PHILLIPS, C. L.; ANDERSON, J. A.; GLOTZER, S. C. Pseudo-random number
generation for Brownian Dynamics and Dissipative Particle Dynamics simulations on GPU
devices. Journal of Computational Physics 230(19): 7191-7201, Agosto 2011.
41.
KAUKONEN, M. et al. Lennard-Jones Parameters for Small Diameter Carbon
Nanotubes and Water for Molecular Mechanics Simulations from van der Waals Density
Functional Calculations. Journal of Computational Chemistry, 2012. 652-659.
42.
MORTENSEN, N. A.; OKKELS, F.; BRUUS, H. Reexamination of Hagen-Poiseuille flow:
Shape dependence of the hydraulic resistance in microchannels. Physical Review E 71,
057301, 2005.
46
43.
WEI, B.; ROGERS, B. J.; WIRTH, M. J. Slip Flow in Colloidal Crystals for Ultraefficient
Chromatography. Journal of the American Chemical Society, Abril 2012.
44.
GÖSCH, M. et al. Hydrodynamic Flow Profiling in Microchannel Structures by Single
Molecule Fluorescence Correlation Spectroscopy. Analytical Chemistry, vol. 72, No 14,
15 Julho 2000. 3260-3265.
45.
FALKENA, W. File Exchange - xml2struct. MatLab Central, 2012. Disponivel em:
<http://www.mathworks.com/matlabcentral/fileexchange/28518-xml2struct>.
Acesso
em: 15 jun. 2013.
46.
HUMPHREY, W.; DALKE, A.; SCHULTEN, K. VMD - Visual Molecular Dynamics. J.
Molec. Graphics, 1996. 33-38.
47
Anexo A – Arquivo de Configuração do Hoomd-blue
A seguir, o arquivo de configuração do Hoomd-blue (método_DPD.hoomd) para que
ele simule o sistema corretamente como especificado. Este exemplo utiliza a força de campo
como 0,0253 uF.
From hoomd_script import *
import random
# -- Parametros para a construcao do paralelepipedo de particulas
a = 1.01 # Porcentagem de espacamento entre particulas
m = 13 # Numero de particulas na largura e altura do paralelepipedo
ld = 3 # Razao de aspecto (L/D)
# criando 12x12x39 particulas moveis e uma caixa de (13x13x39)*1,01 com 1%
de espaco vazio
system = init.create_empty(N=m*m*m*ld, box=(m*a, m*a, m*a*ld),
n_particle_types=3)
# alterando o tipo de particulas para a confeccao do tubo
lo = - m*a / 2.0;
for p in system.particles:
# dividindo as particulas em indices nas direcoes do sistema de
coordenadas
(i, j, k) = (p.tag % m, p.tag/m % m, p.tag/m**2 % (m*ld))
# distribuindo as particulas uniformemente em um paralelepipedo de
# 20x20x60 particulas com espacamento de 1,01 raio entre elas
p.position = (lo + i*a + a/2, lo + j*a + a/2, (lo*ld) + k*a + a/2)
# Ajustando a massa para elevar a densidade
p.mass = 5.151505
# Particulas inicialmente sao do tipo A
p.type = 'A'
# As particulas nas bordas sao do tipo B
if i==0 or i==(m-1) or j==0 or j==(m-1):
p.type = 'B'
# define uma fatia de particulas como sendo do tipo C
# para verificar o sentido do fluxo e se ele esta ocorrendo
if (i>0 and i<(m-1)) and (j>0 and j<(m-1)) and k>=0 and k<=1:
p.type = 'C'
# inicializando as velocidades das particulas para
# tenham distribuicao normal
random.seed(1234);
T = 1.0
px = py = pz = 0.0;
for p in system.particles:
mass = p.mass;
vx = random.gauss(0, T / mass)
vy = random.gauss(0, T / mass)
vz = random.gauss(0, T / mass)
p.velocity = (vx, vy, vz)
# somando o momento total do sistema
px += mass*vx;
48
py += mass*vy;
pz += mass*vz;
# computando a media
px /= m*m*m*ld;
py /= m*m*m*ld;
pz /= m*m*m*ld;
# subtraindo das particulas a media dos momentos para que a soma seja zero
no sistema
for p in system.particles:
mass = p.mass;
v = p.velocity;
# reduzindo a amplitude da aleatoriedade das particulas
p.velocity = ((v[0] - px/mass)*0.1, (v[1] - py/mass)*0.1, (v[2] pz/mass)*0.1);
# incluindo a dinamica do sistema pelo metodo DPD
dpd = pair.dpd(r_cut=1.0, T=1.0, seed=12345)
dpd.pair_coeff.set('A', 'A', A=25.0, gamma = 0.5)
dpd.pair_coeff.set('A', 'B', A=25.0, gamma = 0.5)
dpd.pair_coeff.set('A', 'C', A=25.0, gamma = 0.5)
dpd.pair_coeff.set('B', 'B', A=25.0, gamma = 0.5)
dpd.pair_coeff.set('B', 'C', A=25.0, gamma = 0.5)
dpd.pair_coeff.set('C', 'C', A=25.0, gamma = 0.5)
# separando em grupos os tipos A, B e C
groupA = group.type(name='groupA', type='A')
groupB = group.type(name='groupB', type='B')
groupC = group.type(name='groupC', type='C')
# incluindo a forca constante Fb em todas particulas fluidicas
force.constant(fx=0.0, fy=0, fz=0.0253, group=groupA)
force.constant(fx=0.0, fy=0, fz=0.0253, group=groupC)
# realizando a integracao de 0,002 pelo metodo Velocity-Verlet
integrate.mode_standard(dt=0.002)
integrate.nve(group=groupA)
integrate.nve(group=groupC)
# gerando arquivo de visualizacao online
xml = dump.xml(filename='dump_metodo_DPD.xml', vis=True)
# salvando em um arquivo .dcd as posicoes das particulas a cada 100
periodos
dump.dcd(filename='dump_metodo_DPD.dcd', period=100)
# salvando arquivos de log para verificacao do perfil de velocidade pelo
software MatLab
dump.xml(filename="atm/atoms_xml", period=1000, velocity=True,
position=True, type=True)
# configurando visualizacao online
analyze.imd(port=54321, period=500)
# executando 40.000 iteracoes nesta simulacao
run(4e4)
49
Anexo B – Programa em MatLab
Programa desenvolvido na plataforma MatLab com objetivo de gerar um gráfico do
perfil de velocidades de certas partículas, dado um arquivo XML que contém as posições e
velocidades das partículas em determinada iteração da simulação.
function [filtred_table] = plotPoiseuilleSurf (pos_z, porc_discret, r_cut,
particles_pos, particles_vel)
% Faz um filtro da media da velocidade das particulas para cada
% ponto discretizado desejado gerando a superficie de velocidades
% dos pontos discretizados.
% #### VARIAVEIS DE ENTRADA ####
% pos_z = integer (1) (posicao estcolhida em Z para verificacao do
perfil
%
de velocidades de Poiseuille)
% porc_discret = float [0..1] (1) (porcentagem do tamanho do perfil do
tubo
%
que sera utilizado como subdivisoes
do
%
dominio)
% r_cut = float (1) (distancia maxima de influencia da velocidade das
%
particulas na vizinhanca do ponto em que sera
%
calculado a velocidade do fluido
% particles_pos = float (NumOfParticles x 3) (vetor posicao das
particulas
%
no dominio)
% particles_vel = float (NumOfParticles x 3) (vetor velocidade de cada
%
particula do dominio)
% #### INICIO DO PROGRAMA ####
% agregar as posicoes e as velocidades em uma matriz soh
sistema = [particles_pos,particles_vel];
% filtrar em Z as posicoes e velocidades de acordo com
% a posicao em Z que se deseja saber a superficie de
% Poiseuille e o raio de corte de influencia da
% velocidade das particulas
sistema = sistema(sistema(:,3)>-r_cut & sistema(:,3)<r_cut,:);
% cria tabela final de velocidades
hash_size = floor(1/porc_discret);
filtred_table = zeros(hash_size,hash_size);
% obtem os valores limites
largura_min = ceil(min(sistema(:,1)));
largura_max = floor(max(sistema(:,1)));
altura_min = ceil(min(sistema(:,2)));
altura_max = floor(max(sistema(:,2)));
step_x = (largura_max-largura_min)/hash_size;
step_y = (altura_max-altura_min)/hash_size;
% para cada passo em x e y oriundo da porcentagem de
% discretizacao
for i = 1:hash_size
for j = 1:hash_size
50
% filtrar em X e Y o dominio de acordo com a posicao
% atual de calculo e com o raio de corte
x = largura_min + step_x * (i-1) + step_x/2;
y = altura_min + step_y * (j-1) + step_y/2;
sistema_local = sistema(sistema(:,1)>x - r_cut &
sistema(:,1)<x + r_cut & sistema(:,2)>y-r_cut & sistema(:,2)<y + r_cut,:);
% zera variaveis internas
cont = 0;
soma_vel = 0;
% para cada particula no dominio
for p = 1:length(sistema_local(:,1))
% calculo da distancia da particula atual
dist = norm(sistema_local(p,1:3)-[x,y,pos_z]);
% verificar se a distancia do ponto ateh
% a particula eh igual ou menor que o raio de corte.
if (dist <= r_cut)
% Se verdadeiro, adiciona a velocidade em Z na
% variavel da posicao da particula
soma_vel = soma_vel + sistema_local(p,6);
% Acrescenta 1 no numero de particulas de influencia
cont = cont + 1;
end
end
if (cont>0)
% Finalizado, divide a reducao das velocidades pelo
% numero de particulas que influenciaram
filtred_table(i,j) = soma_vel/cont;
end
end
end
% Plota o grafico da tabela discretizada com o perfil de velocidades
em 3D
steps_width = largura_min+step_x/2:step_x:largura_max-step_x/2;
steps_height = altura_min+step_y/2:step_y:altura_max-step_y/2;
[X1,X2] = meshgrid(steps_width,steps_height);
f = polyFit2D(filtred_table(1:992,:),X1(1:992,:),X2(1:992,:),2,2);
Z = polyVal2D(f,X1(1:992,:),X2(1:992,:),2,2);
figure;
subplot(1,2,1);
surfc(X1(1:992,:),X2(1:992,:),Z,'LineStyle','none');
title('Perfil de velocidade 3D');
xlabel('largura');
ylabel('altura');
zlabel('velocidade do fluido');
hold all;
% Plota a curva de Hagen-Poiseuille na linha horizontal posicionada no
% centro da grid
subplot(1,2,2);
filtred_table_line = filtred_table(:,500);
f = polyfit(steps_width,filtred_table_line',2);
51
Z = polyval(f,steps_width);
plot(steps_width, Z);
title('Perfil de Hagen-Poiseuille');
xlabel('largura');
ylabel('velocidade do fluido');
hold off;
end
52