Análise Matemática do Método SPH via Estudo do Modelo de Ruptura de Barragens Neylan Leal Dias Orientador: Prof. Dr. Messias Meneguete Júnior Co-orientadora: Profa. Dra. Analice Costacurta Brandi Programa: Matemática Aplicada e Computacional Presidente Prudente, Novembro de 2017 UNIVERSIDADE ESTADUAL PAULISTA Faculdade de Ciências e Tecnologia de Presidente Prudente Programa de Pós-Graduação em Matemática Aplicada e Computacional Análise Matemática do Método SPH via Estudo do Modelo de Ruptura de Barragens Neylan Leal Dias Orientador: Prof. Dr. Messias Meneguette Júnior Co-orientadora: Profa. Dra. Analice Costacurta Brandi Dissertação apresentada ao Programa de Pós-Graduação em Matemática Aplicada e Computacional da Faculdade de Ciências e Tecnologia da UNESP para obtenção do tí- tulo de Mestre em Matemática Aplicada e Computacional. Presidente Prudente,17 de Novembro de 2017 D541a Dias, Neylan Leal Análise Matemática do Modelo SPH via Estudo do Modelo de Ruptura de Barragens / Neylan Leal Dias. -- Presidente Prudente, 2017 56 p. Dissertação (mestrado) - Universidade Estadual Paulista (Unesp), Faculdade de Ciências e Tecnologia, Presidente Prudente Orientadora: Messias Meneguette Jr Coorientadora: Analice Costacurta Brandi 1. Smoothed Particle Hydrodynamics. 2. Ruptura de barragem. 3. SPHysics. I. Título. Sistema de geração automática de fichas catalográficas da Unesp. Biblioteca da Faculdade de Ciências e Tecnologia, Presidente Prudente. Dados fornecidos pelo autor(a). Essa ficha não pode ser modificada. Dedico à minha amada esposa Erica Dias aos meus pais Neilza Leal, Edimilson Vales e Ailton Dias e aos meus irmãos Daniela Mayra e Keven Westin Agradecimentos Agradeço acima de tudo a Deus por proporcionar a oportunidade de cursar o mestrado, pela saúde física, mental e autocontrole diante das situações difíceis enfrentadas em cada etapa deste mestrado. Agradeço a minha mãe Neilza Pantoja Leal, pelo seu amor, dedicação, preocupação e apesar de longe, sempre atenciosa em todas etapas da minha vida. Agradeço a minha amada esposa Erica Brito Dias, pela sua dedicação, aceitando o desafio de mudar para uma cidade distante, pelo tempo que fiquei ausente de casa para os estudos do mestrado, e por me aturar em todos os momentos dessa caminhada. Ao Prof. Dr. Messias Meneguette Jr, que além de orientador, um amigo, pela sua paciência, e acreditar no meu potencial, colaborando para meu crescimento pessoal e profissional. Ao Prof. Dr. Carlos Alberto Dutra Fraga Filho, pela sua enorme contribuição para com a construção dos resultados alcançados por essa dissertação, bem como pela sua disponibilidade e paciência. À profa. Dra. Analice Costacurta Brandi, pela sua disponibilidade e ensinamentos durante o mestrado. Ao prof. Dr. Cássio Oishi, cordenador do mestrado, pela sua disposição para resolver alguns problemas ocorridos, e pelos seus ensinamentos na disciplina ministrada. À profa. Dra. Simone de Almeida Delphin, por me mostrar o caminho da ciência, e sempre me acompanhar nas etapas dos meus estudos, e ao meu grande amigo Edson Leal. Aos meus familiares, meus irmãos Keven Westin Leal dos Santos e Daniela Mayra Pantoja Leal, meu pai Ailton Dias Ferreira, minha avó Sebastiana Pantoja Leal, minha sogra Adalgisa, e sogro Ziomar Brito e em especial ao meu Pai de criação Edimilson Vales dos Santos, pelo seu amor, preocupação e auxílios durante meus estudos. Aos meus amigos da PosMAC: Leonardo, Paulo, Rafael, Anderson, Laison, Rhuan, Jorge, Rodrigo, Tânia, Karlla, Maria cecília, Thais, Guilherme, Luciano, em especial às minhas amigas Débora, grande parceira durante as disciplinas, e trabalhos, e à Clícia pelo apoio durante a pesquisa e construção dos resultados; por inúmeras vezes discutimos acerca do tema estudado. O presente trabalho foi realizado com apoio da Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Código de Financiamento 001. � A matemática é o alfabeto com o qual DEUS escreveu o universo�. Pitágoras Resumo O estudo dos métodos sem malha tem sido muito atrativo nos últimos anos, devido a suas capacidades de melhor resolver de�ciências enfrentadas pelos métodos tradicionais que usam malha. O SPH é um dos primeiros métodos sem malha, por isso vários estu- dos têm sido feitos para seu melhor desenvolvimento. Uma das aplicações clássicas em dinâmica dos �uidos é a ruptura de barragens. Uma das problemáticas enfrentadas na aplicação da ruptura de barragem é a recuperação do campo de velocidades e de pressão. O objetivo deste trabalho é analisar as correções aplicadas ao SPH a �m de suavizar as oscilações dos campos de velocidades e do campo de pressões, usando a viscosidade arti�cial e as correções dos �ltros de densidade Shepard e MLS, por meio do software SPHysics. Além disso, foi feita uma comparação dos resultados simulados com resultados de laboratório. As correções aplicadas se mostraram e�cientes na recuperação dos campos propostos. Quando comparados os resultados simulados para o escoamento do �uido com os resultados laboratoriais, foram obtidas boas aproximações, com diferenças percentuais baixas. Palavras-Chave: Smoothed Particle Hydrodynamics, Filtro de densidade Shepard e Moving Least Squares, Ruptura de barragem, SPHysics. Abstract The study of the meshless methods have been very attractive in recent years due to their ability to better address shortcomings faced by traditional mesh based methods. SPH is one of the �rst meshless methods, so several studies have been done for its best development. One of the classic applications in �uid dynamics is the dam break problem, which has been the subject of much research using the SPH method. One of the di�culties faced in the application of dam break is the recovery of velocity and pressure �elds. The objective of this work is to analyze the corrections applied to SPH in order to smooth and better recover the velocity and pressure �elds, using arti�cial viscosity and the corrections of the Shepard and MLS density �lters, through the SPHysics software. In addition, a comparison between simulated and available experimental results is made in order to vali- date the study. The applied corrections were e�cient for recovering the mentioned �elds. The comparisons showed very good agreement, that is, very low percentage di�erences. Keywords: Smoothed Particle Hydrodynamics, Shepard density �lters and Moving Le- ast Squares, Dam Break, SPHysics. Lista de Figuras 2.1 Descrições euleriana e lagrangiana. . . . . . . . . . . . . . . . . . . . . . . 20 2.2 Diferentes tipos de discretização. Fonte:[6] . . . . . . . . . . . . . . . . . . 20 2.3 Partículas dispostas de forma regular (a) e randômica (b). Fonte: [23] . . 21 3.1 Exemplo de uma distribuição de partículas 2-D envolvendo a partícula a. O raio de in�uência é expresso como um múltiplo, k, do comprimento de suavização h. [22] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2 Figura a) mostra todo o domínio computacional e Figura b) células adja- centes às partículas vizinhas [21]. . . . . . . . . . . . . . . . . . . . . . . . 27 3.3 Núcleo Gaussiano e sua derivada dividida pelo fator dimensional αD. [3] . 28 3.4 Núcleo Quadrático e sua derivada dividida pelo fator dimensional αD. [3] . 28 3.5 Núcleo Spline cúbico e sua derivada dividida pelo fator dimensional αD. [3] 29 3.6 Núcleo Quíntico e sua derivada dividida pelo fator dimensional αD. [3] . . 30 4.1 Condição de repulsão. [22] . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2 Evolução de novas células na direção Z, dependendo das partículas do �uido em movimento. [22] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.1 Escoamento de lava, sem correção XSPH (a) e com (b). [18] . . . . . . . . 38 5.2 Estudo da tensão SPH. [20] . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.3 Árvore do diretório. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.4 Esquema da simulação. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 6.1 Con�guração inicial do tanque . . . . . . . . . . . . . . . . . . . . . . . . . 43 6.2 Campo de velocidades do �uidos parte 1. . . . . . . . . . . . . . . . . . . . 44 6.3 Campo de velocidades dos �uidos parte 2. . . . . . . . . . . . . . . . . . . 44 6.4 Pressão com uso da viscosidade arti�cial α = 0.3, com �ltros de densidade de Shepard e MLS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 6.5 Grá�co dos valores da pressão com uso da viscosidade arti�cial α = 0.3 e com �ltros de densidade de Shepard e MLS. . . . . . . . . . . . . . . . . . 46 6.6 Geometria do tanque para validação experimental. [4] . . . . . . . . . . . . 47 6.7 Geometria 1 (a) e geometria 2 (b) para simulações. . . . . . . . . . . . . . 47 6.8 Comparação com modelo experimental, para α = 0.2 . . . . . . . . . . . . 48 6.9 Comparação com modelo experimental, para α = 0.3 . . . . . . . . . . . . 49 6.10 Comparação com modelo experimental, segunda geometria, para α = 0.2 . 49 6.11 Comparação com modelo experimental, segunda geometria, para α = 0.3 . 50 6.12 Dimensões do tanque 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 6.13 Comparação da interface do �uido simulada com experimento laboratorial, para os tempos 0.01s e 0.16s . . . . . . . . . . . . . . . . . . . . . . . . . 51 6.14 Comparação da interface do �uido simulada com experimento laboratorial, para os tempos 0.27s e 0.37s . . . . . . . . . . . . . . . . . . . . . . . . . 51 9 LISTA DE FIGURAS 10 6.15 Comparação da interface do �uido simulada com experimento laboratorial, para o tempo 0.45s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 Lista de Tabelas 1.1 Rompimento de barragens no Brasil. . . . . . . . . . . . . . . . . . . . . . 16 6.1 Valores médios da pressão de acordo com o tempo, referentes as análises da �gura (6.5) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 6.2 Resultados das �guras (6.8) e (6.9), primeira geometria. . . . . . . . . . . . 48 6.3 Resultados das �guras (6.10) e (6.11), segunda geometria. . . . . . . . . . . 48 6.4 Resultados das �guras (6.13) e (6.14), segunda geometria. . . . . . . . . . . 52 Sumário Resumo 5 Abstract 7 Lista de Figuras 8 Lista de Tabelas 10 Capítulos 1 Introdução 15 1.1 Aplicações SPH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.2 Ruptura de Barragens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.3 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2 Matemática do Fluido 19 2.1 Descrição do Fluido . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2 Modelagem Computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3 Fundamentos do SPH 23 3.1 Interpolação Padrão SPH . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1.1 Erros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Núcleo de Suavização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2.1 Busca por partículas vizinhas . . . . . . . . . . . . . . . . . . . . . 26 3.3 Tipos de Núcleos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4 Modelagem Matemática 31 4.1 Equação de Quantidade de Movimento . . . . . . . . . . . . . . . . . . . . 31 4.2 Equação de Estado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.3 Conservação de Massa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.4 Conservação de Energia . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.5 Tratamento de Fronteira . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.5.1 Fronteira Sólida . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.5.2 Fronteira Livre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 5 Correções Aplicadas ao SPH 35 5.1 Viscosidade Arti�cial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.2 Filtro de Densidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 5.2.1 Filtro Shepard - Ordem Zero . . . . . . . . . . . . . . . . . . . . . . 36 5.2.2 Mínimos Quadrados Móveis (MLS) - Primeira Ordem . . . . . . . 36 5.3 Renormalização do Núcleo . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.3.1 Correção do Gradiente do Núcleo . . . . . . . . . . . . . . . . . . . 37 5.4 Movimento das Partículas . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.5 Correção da Tensão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.6 Atualização Temporal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.7 Simulador . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 6 Resultados 43 6.1 Campo de velocidades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 6.2 Campo de pressões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 6.3 Frente de onda: Validação experimental, tanque 1. . . . . . . . . . . . . . . 46 6.4 Frente de onda: Validação experimental, tanque 2. . . . . . . . . . . . . . . 49 7 Conclusão 53 Referências 53 Apêndices Capítulo 1 Introdução Desde o seu nascimento em 1977, o Smoothed Particle Hydrodynamics (SPH) tem sido melhorado, atingindo resultados promissores. Isso se deve a um intenso esforço teó- rico pela comunidade trabalhando nesta area e também a uma evolução na capacidade computacional, principalmente em computação paralela. Por ser um método sem malha e relativamente recente já é considerado uma ferramenta promissora em aplicações industri- ais e ambientais, pois pode ser usado com certa facilidade para problemas com geometrias complexas e fronteiras livres. Nos últimos anos, alguns códigos SPH tornaram-se uma parte inerente do arsenal nu- mérico de laboratórios industriais de pesquisa e desenvolvimento e instituições acadêmicas. Existem várias razões que explicam este sucesso [22]: (1) a crescente necessidade da in- dústria e pesquisa em utilizar ferramentas apropriadas para solução numérica de modelos hidrodinâmicos complexos; (2) avanços recentes na teoria SPH que resolveram uma série de problemas deste método; e (3) o surgimento das unidades de processamento grá�co (GPUs) que permitem aos códigos SPH a simulação de �uxos tridimensionais complexos em escalas reais, mantendo tempos computacionais razoáveis. 1.1 Aplicações SPH O SPH foi originalmente desenvolvido para a modelagem astrofísica, onde a densidade do �uido pode variar no espaço e com o tempo em algumas ordens de grandezas e na ausência de fronteiras físicas. A capacidade do SPH lidar com as questões de densidade variável e �uxos diversos, decorre do fato de sua natureza ser totalmente lagrangiana. A idéia do método é acom- panhar o �uido por meio das partículas em movimento, sem qualquer restrição de malha. Em meados da década de 1990 Monaghan [12] fez uma primeira tentativa de usar SPH para �uidos simples, ou seja, aplicou-se à hidráulica de superfície livre encontrada na vida cotidiana. Para isto, Monaghan apresentou duas idéias básicas, mas úteis, para resolver os problemas: aproximar densidade constante e impor limites rígidos. A idéia de impor uma densidade constante equivale a uma condição quase incompressível, ou seja, uso da equação de estado com velocidade de som su�cientemente grande que não impedisse simu- lações. Esta é a variante conhecida atualmente como WCSPH (SPH quase compressível). Surpreendentemente, o tratamento da superfície livre não exigiu cuidado especial nessa tentativa inicial. Este primeiro sucesso revelou-se tão impressionante que a comunidade de modeladores do SPH na hidráulica cresceu rapidamente desde meados da década de 1990, propondo novas idéias e melhorias para o modelo simples original de Monaghan. 15 1. Introdução 16 O sucesso do SPH na hidráulica é devido à sua extraordinária capacidade de simular uma variedade de �uxos complexos, envolvendo um amplo aspecto de processos físicos. As principais características e mais interessantes [22] permitem lidar com : • superfícies livres com deformações acentuadas e rápida movimentação sem qualquer restrição em relação às suas topologias; • �uxos altamente não-lineares, dominados por inércia e processos de impacto; • �uxos multi-�uidos e multi-físicos; • interações �uido / estrutura, incluindo o movimento rígido do corpo em um acopla- mento �uido e �uido / elasticidade; • facilidade de programação (pelo menos para o básico), seja em 2-D ou 3-D. Uma das caracteristicas atraentes do SPH é seu caráter variacional, que mantem har- monia com os fundamentos da mecânica lagrangiana e hamiltoniana (para fenômenos não dissipativos). Desde seu uso na hidráulica, o SPH foi aplicado a diversos problemas na pesquisa de escoamentos e engenharia hidráulica. O sucesso deste método em vários cam- pos deve-se à relativa facilidade com que as simulações SPH foram capazes de produzir resultados para casos não-lineares complicados e muitas vezes fenômenos multi-fásicos. Com pouca modi�cação da metodologia básica, o SPH conseguiu gerar resultados em estreita concordância com soluções de referência/dados em testes de validação, sem algo- ritmos altamente so�sticados exigidos em esquemas baseados em malha. A maioria dos campos acima mencionados são considerados muito difíceis (para não dizer impossíveis) para outros métodos numéricos. No modelo de escoamento com superfície livre, o SPH tem rivalizado com o método dos volumes �nitos e de outras abordagens dedicadas a esse tipo especial de escoamento. Além disso, o progresso recente em sanar aspectos numéricos antes de�citários do SPH, aumentam a credibilidade no uso do método e faz crescer a comunidade que se dedica a esta abordagem. 1.2 Ruptura de Barragens O rompimento de barragens é uma modalidade de desastres consideravelmente reinci- dente na história da humanidade. Dois são os principais fatores que podem ser apontados como causa primária desse evento: o advento de um fenômeno natural intenso respon- sável por abalar a estrutura da barragem ou o mau planejamento dessa estrutura que independentemente de fatores externos entra em colapso em razão dos erros de cálculos dos engenheiros [1]. No Brasil há ocorrência de rompimento de represas com dejetos minerais e também em barragens de água nas últimas décadas, com vítimas, como observado na Tabela (1.1). Tabela 1.1: Rompimento de barragens no Brasil. Local Ano Danos causados Itabirinto 1986 7 óbitos Nova Lima 2001 5 óbitos Miraí 2007 4000 pessoas desabrigadas. Mariana 2015 19 óbitos, 8 desaparecidos, 600 desabrigados. O problema da quebra de barragens foi introduzido pela primeira vez por Stoker (1957). Este problema consiste em ter um espaço fechado cheio de água; a barreira para um lado é então removida e a água pode �uir livremente para dentro do vazio. Para essa simulação é utilizado um código 2D-SPH (2D-SPHysics). A ruptura da barragem é simulada com diferentes combinações de opções de compilação que modi�cam 1. Introdução 17 tratamento numérico e teóricos da modelagem. Assim, os diversos resultados do método SPH são comparados entre si e também contra os dados experimentais existentes, para se encontrar uma combinação adequada que forneça a melhor precisão para a simulação. 1.3 Objetivos O objetivo deste trabalho, é apresentar uma análise matemática do método SPH atra- vés do estudo do modelo de ruptura de barragens. Este trabalho está dividido da seguinte forma: o primeiro capítulo traz a introdução, onde é exposto de forma resumida o pro- blema e aspectos gerais da evolução do método SPH e sua utilização; o segundo capítulo versa sobre alguns conceitos básicos sobre �uidos, que são importantes para a compre- ensão da aplicação estudada neste trabalho; o terceiro trata dos aspectos fundamentais sobre o método SPH, através de sua formulação matemática, erro e apresentação dos nú- cleos; o quarto capítulo aborda as equações de conservação e governantes do escoamento juntamente com suas discretizações SPH; o capítulo cinco descreve as correções que são aplicadas ao SPH, de tal forma que o método possa modelar o �uido corretamente; o ca- pítulo seis traz os resultados obtidos nas simulações e, por último, o capítulo sete �naliza o trabalho com as conclusões. Capítulo 2 Matemática do Fluido Pode-se descrever várias características a respeito do �uido, dentre elas, apresenta-se a seguir as principais. 2.1 Descrição do Fluido É interessante descrever sumarizadamente as concepções euleriana e lagrangiana de um escoamento. Assim, as quantidades físicas associadas aos elementos de �uido variam em um meio contínuo e podem ser descritas como dois tipos [19]: Euleriana Ao invés de acompanhar o movimento ao longo do tempo do escoamento, considera-se um conjunto de posições �xas no espaço como se fossem partículas e, em seguida calcula- se as variações das quantidades físicas do �uido sempre nessas posições. Lagrangiana O escoamento é representado por uma coleção de partículas de �uido e cada uma delas se locomove no tempo, muda de posição no espaço, e então sobre cada partícula são calculadas as novas quantidades físicas do �uido. Observe que elas poderão estar em posições espaciais diferentes em cada tempo; por exemplo, como na �gura (2.1). A �gura (2.1) procura dar ideia dos dois tipos de abordagens comumente usados na simulação numérica. O método SPH utiliza a abordagem lagrangiana. Compressibilidade Um �uido que tem um grande variação em sua densidade é de�nido como compres- sível, já quando não há variação (�xo) é chamado de incompressível. Considera-se um �uido ligeiramente compressível, ou quase-compressível, por exemplo a água, quando há variações de densidades muito baixas. Neste trabalho iremos tratar deste último caso. 19 2. Matemática do Fluido 20 Figura 2.1: Descrições euleriana e lagrangiana. 2.2 Modelagem Computacional O uso de computadores para a simulação computacional de modelos matemáticos é fundamental para a evolução da engenharia, ciência e tecnologia. Estas simulações envolvem experimentos em laboratórios e computações numéricas baseadas em modelos matemáticos. Para efetuar a simulação do modelo matemático é preciso discretizar o domínio por meio de malha estruturada ou não, ou por meio de partículas, ou seja, sem malha. Na �gura (2.2) podemos ver diferentes casos de discretizações. Figura 2.2: Diferentes tipos de discretização. Fonte:[6] Modelos com malha: Malhas computacionais são constituídas por linhas e pontos. Os pontos são conside- rados onde essas linhas se interceptam e servem de posição para o cálculo de propriedades físicas conforme o modelo matemático. Essas malhas podem ser malhas, grades, volumes, células ou elementos, de acordo com o método utilizado. Muitas utilizam pontos igual- mente espaçados, como é o caso de diferenças �nitas [5], outras necessitam de adaptações de acordo com o domínio adotado. 2. Matemática do Fluido 21 Sem Malha: Métodos sem malha (ou meshfree) usam um conjunto de partículas espalhadas dentro do domínio do problema. Pensando cada partícula como um nó, existe iteração entre cada nó com todos os seus vizinhos, assim os nós carregam propriedades como velocidade, energia, massa, etc. As partículas podem estar distribuídas de forma (�xa) regular ou randômica [23]; o mais comum é que as partículas se movimentem (mudando de posição espacial) com o tempo. A �gura 2.2 ilustra as situações: com malha e sem malha. Já a �gura (2.3) exempli�ca as distribuições uniforme e randômica. Figura 2.3: Partículas dispostas de forma regular (a) e randômica (b). Fonte: [23] Capítulo 3 Fundamentos do SPH O Smoothed Particle Hydrodynamics (SPH) é um método lagrangiano que não requer malha ou grade. O material é representado por elementos de massa, chamados de partí- culas. A simulação avança no tempo calculando a nova posição e respectivas propriedades de cada partícula, por meio de um esquema de integração. 3.1 Interpolação Padrão SPH Pode-se conceituar, do ponto de vista computacional, que o �uido é representado po uma porção de partículas evoluindo com a velocidade do escoamento. A base da formu- lação SPH é uma interpolação integral que utiliza a função δ de Dirac, introduzida pelo físico teórico Paul Dirac em 1930. Trata-se de uma função generalizada ou distribuição geralmente representado por δ tal que δ(r− r′) = { 1, se r = r′, 0, se r 6= r′. Na modelagem computacional, o �uido é apresentado como uma parcela de partículas que evoluem com o movimento do escoamento. Cada partícula representa um ponto de interpolação no qual todas as propriedades do �uido são conhecidas. Assim, considerando uma função A de um campo escalar, em um domínio Ω ⊆ R, com posições r,r' e pela de�nição da distribuição do delta de Dirac δ, é possível escrever A(r) = ∫ Ω A(r′)δ(r− r′)dr′. (3.1) Devido a função delta de Dirac ser generalizada, necessita de propriedades de continuidade e diferenciabilidade, logo não pode ser empregada para modelos discretizados numerica- mente. A partir disso, aproxima-se δ(r− r′) através de uma função suave W (r− r′, h) [8] parametrizada por um comprimento de suavização h, que de�ne a área de in�uência da função núcleo W . É exigido então∫ Ω W (r− r′, h)dr′ = 1, lim h→0 W (r− r′, h) = δ, onde o limite é pretendido no sentido das distribuições, onde r é a posição do vetor; W é a função de ponderação ou núcleo, h é chamado de comprimento de suavização e controla 23 3. Fundamentos do SPH 24 o domínio de in�uência Ω. A�m de aproximar o valor da função A(r) coloca-se A(r) = ∫ Ω A(r′)W (r− r′, h)dr′. (3.2) A formulação SPH da integral na equação (3.2) em sua forma discreta é representada por um somatório sobre todas as partículas no domínio de in�uência da partícula a, ou seja, com as partículas b vizinhas mais próximas. O domínio contendo as partículas vizinhas é chamado de domínio do suporte compacto de W . Em termos físicos o domínio Ω representa o corpo do �uido, que será discretizado com um conjunto �nito de partículas, cada uma com massa ma e densidade ρa. O volume in�nitesimal dr′ na integral (3.2) é substituído pelo volume da partícula bth que será identi�cado por ∇Vb e dado pela expressão: ∇Vb = mb ρb . (3.3) Na notação discreta consideramos Wab = W (ra − rb, h). Assim a função A pode ser aproximada usando o valor das partículas próximas por um somatório: A(ra) ≈ ∑ b W (ra − rb, h)A(rb) mb ρb . (3.4) Esta é a formulação básica do SPH, que calcula o valor da função por interpolação integral. Formalmente, seja uma função de densidade ρ : Ω→ R e um conjunto de pontos {ra} ⊂ Ω. Assim a aproximação em notação discreta, leva à seguinte aproximação da função em uma partícula a. A(ra) = ∑ b mb Ab ρb Wab (3.5) onde o somatório considera todas as partículas dentro da região do suporte compacto da função núcleo centrada em a. 3.1.1 Erros Baseado na equação (3.4) podemos encontrar a fórmula diretamente para um gradi- ente: ∇A(ra) = ∑ b mb ρb A(rb)∇bWab (3.6) onde o termo ∇bWab indica o gradiente do núcleo de suavização em relação as coordenadas da partícula a: ∇aWab = ( ∂ ∂xa i + ∂ ∂ya j + ∂ ∂za k ) Wab, (3.7) onde i, j e k são vetores unitários nas suas coordenadas direcionais. No entanto usando a equação (3.6) para uma função A pode acontecer uma resposta sem signi�cado físico para o gradiente, levando a erros signi�cativos [11]. Por esta razão, é sempre sugerido na literatura o uso da seguinte identidade matemática para calcular o gradiente SPH: ρ∇A = ∇(ρA)− A∇ρ. (3.8) Essa identidade leva a seguinte formulação para o gradiente SPH: ∇Ab = ∑ b mb ρb (Ab − Aa)∇aWab. (3.9) 3. Fundamentos do SPH 25 A modelagem das forças entre as partículas tanto na primeira como na segunda fórmula (equações (3.6) e (3.9)) dão uma reação igual e oposta [11]. Para criar um gradiente mais simétrico, a seguinte identidade pode ser usada ∇ ρ = ∇( A ρ )− A ρ2 ∇ρ. (3.10) Essa identidade leva à seguinte formulação para o gradiente SPH: ∇Aa = ∑ b mb ( Ab ρ2 b + Aa ρ2 b ) ∇aWab. (3.11) O erro da interpolação integral pode ser encontrado usando uma função de expansão da série de Taylor A(r) em torno de r [11], cujos dois primeiros termos são: A(r′) = A(r) + (r′ − r).∇A(r) +O ( |r′ − r|2 ) . (3.12) Para o último termo, o erro de truncamento de segunda ordem, a ordem da distância entre as posições r e r' é geralmente semelhante à ordem do comprimento de suavização h. Se usarmos a interpolação SPH na equação 3.12 podemos obter: A(r) = A(r) ∫ Ω W (r− r′, h)dΩ +∇A(r) ∫ Ω (r′ − r)W (r− r′, h)dΩ +O(h2). (3.13) O núcleo de suavização é uma função uniforme em relação a r e, portanto, (r′− r)W (r− r′, h) é uma função impar [8] tornando o segundo termo igual a zero. Além disso, devido às propriedades do núcleo, que são descritas na próxima seção, a primeira integral é igual à unidade, levando-nos a: A(r) = A(r) +O(h2). (3.14) Comprovamos assim que a interpolação SPH é pelo menos de segunda ordem de precisão. Este resultado pressupõe que as integrais podem ser aplicadas a todo o volume especi�cado pelo comprimento de suavização. Perto da fronteira, isso não se aplica, levando a uma precisão reduzida e possíveis instabilidades numéricas. 3.2 Núcleo de Suavização O desempenho de um modelo SPH depende da escolha das funções de suavização. Esta deve tender para a função delta à medida que o comprimento de suavização se aproxima de zero: lim h−→0 W (r− r′, h) = δ(r− r′) (3.15) Os núcleos de suavização são estendidos apenas para partículas em uma vizinhança de a, por isso estes precisam possuir suporte compacto como na �gura(3.1). O núcleo precisa obedecer as seguintes propriedades: 1. su�cientemente suave: W ∈ C > k, k > 1 2. normalizado: ∫ R W (r− r′, h)dr′ = 1 3. Fundamentos do SPH 26 Figura 3.1: Exemplo de uma distribuição de partículas 2-D envolvendo a partícula a. O raio de in�uência é expresso como um múltiplo, k, do comprimento de suavização h. [22] 3. suporte compacto: W (r− r′, h) = 0, quando: |r− r′| = ~rab > kh 4. positivo: W (r− r′, h) ≥ 0 5. decrescente: W (r− r′, h) < W (r− r′, h), se: |r− r′| > |ra − rk| 6. simétrico: W (r− r′, h) = W (r′ − r, h) = W (~r) 7. convergente: lim h→0 W (r− r′, h) = δ(r− r′). A terceira condição é fundamental pois permite a aproximação ser gerada localmente. O domínio local corresponde à zona de in�uência do núcleo. Devido ao suporte compacto, o esforço computacional é reduzido ao usar apenas as partículas vizinhas. 3.2.1 Busca por partículas vizinhas Devido à movimentação das partículas, é preciso manter o controle sobre as partículas. A�m de explorar as vantagens do suporte compacto do SPH é usado um algoritmo de busca de vizinhas, baseado na subdivisão do domínio em células quadradas, de lado 2h. Assim é possível restringir a pesquisa das partículas vizinhas apenas às células adjacentes como mostrado na �gura (3.2). Existem várias funções que obedecem as propriedades mostradas. Para garantir que a integral do núcleo seja igual à unidade é usado um fator de normalização αD, descrito na 3. Fundamentos do SPH 27 Figura 3.2: Figura a) mostra todo o domínio computacional e Figura b) células adjacentes às partículas vizinhas [21]. próxima subseção, que varia de acordo com as dimensões do domínio adotado. O núcleo é então escrito como função de uma variável não dimensional que é a distância da partícula ao longo do comprimento de suavização. Por exemplo, a variavel em uma dimensão �ca: q = |r− r′| h . (3.16) 3.3 Tipos de Núcleos A seguir iremos descrever os núcleos que serão utilizados nas simulações. Várias apro- ximações de núcleo são descritas em [8] e [11]. Gaussiano O núcleo Gaussiano é uma função suave, mesmo para derivadas de ordem mais altas [8]. Por ser uma série in�nita, possui um alto custo computacional e não possui suporte compacto. É representado por, W (r, h) = αDe −q2 , (3.17) onde αD =  1/(h √ π), para 1D, 1/(πh2), para 2D, 1/( √ pih3, para 3D. (3.18) Quadrático Esta função foi utilizada para simular o problema de impacto de alta velocidade [3], sendo que foi então a�rmado que esta função evita o agrupamento de partículas. Possui um custo computacionalmente baixo. Não tem nenhum extremo em seu gradiente, mas tem uma precisão reduzida por ser um polinômio do segundo grau. W (r, h) = α [ 3 16 q2 − 3 4 q + 3 4 ] (3.19) 3. Fundamentos do SPH 28 Figura 3.3: Núcleo Gaussiano e sua derivada dividida pelo fator dimensional αD. [3] onde αD =  3/4h, para 1D, 2/(πh2), para 2D, 5/(4πh3), para 3D. (3.20) Figura 3.4: Núcleo Quadrático e sua derivada dividida pelo fator dimensional αD. [3] Spline cúbico É o mais utilizado, este núcleo é uma aproximação do núcleo gaussiano, tornando-se computacionalmente barato [14], e possui seu suporte compacto. A equação do Núcleo Spline cúbica é da forma W (r, h) = αD { 1− 3 2 q2 + 3 4 q3, se 0 ≤ q ≤ 1, 1 4 (2− q)3, se 1 ≤ q ≤ 2 (3.21) 3. Fundamentos do SPH 29 onde αD =  2/3h, para 1D, 10/(7πh2), para 2D, 1/(πh3), para 3D. (3.22) Figura 3.5: Núcleo Spline cúbico e sua derivada dividida pelo fator dimensional αD. [3] Quíntico É um polinômio de alta ordem e pode capturar efeitos de ordem superior com precisão aprimorada. Essa função é a que apresenta melhor custo benefício, entre precisão e custo computacional. Em geral, quanto maior a ordem dos núcleos, maior a precisão para a aproximação SPH [3]. W (r, h) = αD(1− q 2 )4(2q + 1) 0 ≤ q ≤ 2 (3.23) onde αD =  3/4h, para 1D, 7/(4πh2), para 2D, 7/(8πh3), para 3D. (3.24) 3. Fundamentos do SPH 30 Figura 3.6: Núcleo Quíntico e sua derivada dividida pelo fator dimensional αD. [3] Capítulo 4 Modelagem Matemática Neste capítulo são apresentadas as discretizações das equações de modelagem de �ui- dos, usando SPH. Considera-se que o �uido tem uma baixa compressibilidade, e portanto, se pode usar uma equação de estado para a relação entre pressão e densidade. 4.1 Equação de Quantidade de Movimento A equação de conservação do momentum é dada por dv dt = −1 ρ ∇P + g + Γ, (4.1) onde: Γ são os termos difusivos relacionados as forças de viscosidade, g=(0, 0,−9.81)ms−2 é a aceleração gravitacional, ρ é a massa especí�ca, v é a velocidade, P é a pressão e ∇P é a força devido a pressão, por unidade de volume. A equação (4.1) é descrita para um �uido Newtoniano. A equação do momento ge- ralmente é tratada em sua forma inviscida, assim para se ter em conta o termo viscoso, é utilizado a viscosidade arti�cial, que será mostrada no capítulo 6. 4.2 Equação de Estado O �uido compressível é aproximado a um �uido incompressível por meio de um �uido quase-compressível, onde a pressão é calculada através de uma equação de estado, o que torna o processo mais rápido do que o tradicional que utiliza a equação de Poisson. Outra limitação da compressibilidade é imposta pelo fato de a velocidade do som ser cerca de dez vezes maior do que a velocidade máxima do �uido, mantendo assim variações de densidade no máximo 1%. Assim a pressão é calculada através da equação de estado conhecida como equação de Tait, P = B(( ρ ρ0 )γ − 1), (4.2) onde: B = c2ρ0/γ, é o termo relacionado às �utuações das massas especí�cas do �uido, 31 4. Modelagem Matemática 32 ρ0 é a massa especí�ca de repouso do �uido e γ = 7. c é a magnitude da velocidade de propagação do som no �uido. 4.3 Conservação de Massa Como o �uido aqui tratado é a água, no SPH, é considerado como fracamente com- pressível (ou quase-compressível). Assim, não é possível garantir que a massa especí�ca seja constante. Considerando a equação de conservação de massa dρ dt = −ρ~∇.~v. (4.3) No formalismo padrão do SPH o �uido é tratado como compressível, o que permite o uso de uma equação de estado para determinar a pressão do �uido, em vez de resolver outra equação diferencial. No entanto, a compressibilidade é ajustada para dimunuir a velocidade do som para que o tempo do modelo (com base na velocidade do som) seja razoável. Da equação (4.3) na notação SPH, calcula-se as mudanças na massa especí�ca do �uido por meio de: dρa dt = ∑ b mbvab∇aWab (4.4) em vez de usar uma soma ponderada de termos de massa, pois é sabido resultar em uma diminuição de densidade arti�cial em relação a interfaces de �uidos [11]. Pode-se também discretizar a densidade como uma função, �cando na forma como na equação (4.5), que é utilizada em problemas onde não há mudança considerável da massa especí�ca; considera- se o �uido quase-compressível [14]. < ρi >= N∑ b=1 mbWab. (4.5) 4.4 Conservação de Energia Durante a simulação, calcula-se a energia cinética, potencial e térmica. A energia térmica associada a cada partícula usando viscosidade arti�cial é dada por de dt = 1 ρ − P · ∇v + εv +∇ · q + qH (4.6) que na notação SPH é calculada pela expressão dada por [11] dea dt = 1 2 ∑ b ( Pa ρ2 a + Pb ρ2 b + ∏ ab ) vab∇aWab (4.7) onde ∏ ab é o termo de viscosidade arti�cial. A energia total do sistema é calculada como a soma da energia cinética e potencial. 4.5 Tratamento de Fronteira A imposição de fronteiras sólidas impermeáveis é essencial para a aplicação de ruptura de barragens. Dentre as existentes na literatura, escolheu-se a de funções repulsivas. Pode também ocorrer o desprendimento de partículas do �uido (fronteira livre) que no SPHysics é identi�cado e tratado adequadamente, como observado na �gura(4.2). 4. Modelagem Matemática 33 4.5.1 Fronteira Sólida As partículas do �uido sofrem uma força repulsiva cuja magnitude está em função direta com a distância entre a partícula e a fronteira, como mostrado na �gura (4.1). Em [12] é utilizada uma função de repulsão para de�nir uma força que mantem as partículas dentro do domínio computacional. No caso de uma partícula que se move paralelamente à fronteira o resultado não é satisfatório [15]. Porém mais tarde, em [15] esta técnica foi modi�cada para se adequar às partículas que se movem paralelamente a uma parede. Figura 4.1: Condição de repulsão. [22] 4.5.2 Fronteira Livre No SPH, as partículas do �uido podem se desprender de diferentes formas. Essas partículas devem ser identi�cadas e tratadas para evitar efeitos espúrios. O tratamento dessas partículas é chamado, no SPHysics, de "Veri�cação dos Limites"[3]. Quando uma partícula do �uido supera o limite superior na direção Z vertical, o domí- nio computacional é estendido e novas células são criadas, como observado na �gura (4.2). As partículas do �uido podem então ocupar essas células de acordo com sua posição. As- sim, o número de células na vertical é modi�cado dinâmicamente dependendo da posição da "maior"partícula do �uido. O caso inverso ocorre quando as partículas descem, como mostrado na �gura (4.2). É interessante observar que o SPHysics também faz esse mesmo controle nas direções X e Y, para aumentar a e�ciência computacional [3]. 4. Modelagem Matemática 34 Figura 4.2: Evolução de novas células na direção Z, dependendo das partículas do �uido em movimento. [22] Capítulo 5 Correções Aplicadas ao SPH Embora a cinemática das simulações de SPH seja geralmente realista, o campo de pressão das partículas pode exibir grandes oscilações, que podem ser particularmente importantes nas áreas próximas de superfícies livres. Um dos métodos mais simples e computacionalmente menos onerosos para suavizar as oscilações do campo pressão, mesmo no caso do uso da equação de Tait, é a utilização de um termo de viscosidade arti�cial, conforme a equação (5.1). Outra técnica neste mesmo sentido é o �ltro de Shepard, também abordado neste capitulo. 5.1 Viscosidade Arti�cial Em problemas que envolvem ondas de choque, ocorre a transformação da energia cinética em calor, é representada como uma forma de dissipação viscosa. Seu uso tem por �nalidade diminuir as instabilidades numéricas e a interpenetração entre partículas. A equação (4.1) pode ser escrita da forma: dva dt = ∑ b mb ( Pa ρ2 a + Pb ρ2 b + ∏ ab ) ∇aWab + g (5.1) O termo gradiente de pressão em forma simétrica é expresso em notação SPH como −1 ρ ∇P = − ∑ b mb ( Pa ρ2 a + Pb ρ2 b ) ∇aWab (5.2) onde P é a pressão e ρ a densidade. O termo ∏ ab é o de viscosidade arti�cial: ∏ ab=  −αcabµab ρab se vab.rab < 0, 0, c.c.. com µab = hvab.rab r2 ab + η2 (5.3) ρab = 1 2 (ρa + ρb), cab = 1 2 (ca + cb) e η2 = 0.01h2, 35 5. Correções Aplicadas ao SPH 36 onde: α um parâmetro livre que pode ser alterado de acordo com cada problema. ca e cb são as velocidades do som nas partículas �xas e vizinhas respectivamente, ha e hb são os raios de suporte das partículas �xas e vizinhas, respectivamente, ϕ2 é um fator que evita diferenças numéricas quando duas partículas se aproximam. O termo relacionado à viscosidade arti�cial é adicionado aos termos de pressão nas aproximações de SPH para as equações de equilíbrio de momento e conservação de energia [7]. 5.2 Filtro de Densidade As partículas no SPH podem sofrer grandes variações de pressão, principalmente nas proximidades das fronteiras e da superfície livre. Para diminuir o efeito da variação, foram introduzidas algumas correções. Um dos problemas é que o número de partículas varia nessas regiões mais afetadas. Assim, uma das correções busca adequar o valor de cada massa especí�ca que ocorre no suporte compacto do núcleo, isso corresponde a uma "�ltra- gem" dos valores em questão. Assim, corrigindo por meio de um �ltro de densidade, nas massas especí�cas das partículas vizinhas, consegue-se suavizar computacionalmente as �utuações de pressão. Esse procedimento pode ser repetido periodicamente "reiniciando" as propriedades em intervalos de�nidos de iterações. 5.2.1 Filtro Shepard - Ordem Zero O �ltro Shepard é uma correção rápida e simples para o campo de densidade. O seguinte procedimento é aplicado a cada m especi�cado entre 20 e 50 passos do tempo ρnewa = ∑ b ρjŴab mb ρb = ∑ b mbŴab (5.4) onde o núcleo foi corrigido usando uma correção de ordem zero, Ŵab = Wab∑ bWab mb ρb . (5.5) Esta técnica possui um baixo custo computacional. 5.2.2 Mínimos Quadrados Móveis (MLS) - Primeira Ordem A abordagemMoving Least Squares (MLS) foi desenvolvida por Dilts (1999) e aplicada por Colagrossi e Landrini (2003). O método é uma correção de primeira ordem para que a variação linear do campo de densidade possa ser reproduzida exatamente: ρnewa = ∑ ρbW MLS ab mb ρb = ∑ b mbW MLS ab (5.6) 5. Correções Aplicadas ao SPH 37 O núcleo corrigido é calculado da seguinte forma: WMLS ab = β(xa.(xa − xb))Wab (5.7) De modo que, por exemplo, em 2-D WMLS ab = [β0(xa) + β1x(xa)(xa − xb) + β1z(xa)(za − zb)]Wab, (5.8) onde o vetor de correção β é dado por β(xa) =  β0 β1x β1z  = A−1 =  1 0 0  (5.9) A = ∑ b Wab mb ρb (5.10)  =  1 (xa − xb) (za − zb) (xa − xb) (xa − xb)2 (za − zb)(xa − xb) (za − zb) (xa − xb)(za − zb) (za − zb)2  . (5.11) Semelhante ao �ltro Shepard, a correção de densidade deve ser aplicada a cada passo do tempo. Nota-se que na aplicação do interpolador MLS, se faz necessária a solução de um sistema linear, em que a matriz dos coe�cientes  pode ser singular, devido à falta de partículas vizinhas. Assim, faz-se uma veri�cação da possibilidade de redução da matriz, para uma submatriz com determinante não nulo, como uma propriedade de envelopamento da matriz dos coe�cientes, uma das características do interpolador MLS. 5.3 Renormalização do Núcleo A correção periódica da função núcleo W é necessária nos cálculos hidráulicos SPH, onde o domínio �nito e a superfície livre são muitas vezes parte do domínio computacional. As partículas perto dos limites ou na superfície livre possuem uma função de alisamento do núcleo truncada devido à ausência de partículas vizinhas. A condição de consistência de W falha. Corrigindo oportunamente a função W ou o seu gradiente resolve-se este problema. 5.3.1 Correção do Gradiente do Núcleo Garante que o gradiente do campo de velocidades seja avaliado corretamente, modi�- cando o gradiente do núcleo através da introdução de uma matriz de correção L: ∇̃Wb(ra) = La∇Wb(ra) (5.12) O gradiente da velocidade é consequentemente estimado como ∇va = N∑ b=1 mb ρb (vb − va)⊗ ∇̃Wb(ra) (5.13) = N∑ b=1 mb ρb (vb − va)⊗ La∇Wb(ra). (5.14) 5. Correções Aplicadas ao SPH 38 O gradiente do Núcleo corrigido deve satisfazer a condição N∑ b=1 mb ρb (rb − ra)⊗∇Wb(ra) = I. (5.15) Assim N∑ b=1 mb ρb (rb − ra)⊗ ∇̃Wb(ra) (5.16) = ( N∑ b=1 mb ρb (rb − ra)⊗∇Wb(ra) ) LTa = I (5.17) da qual L é avaliado como La = ( N∑ b=1 mb ρb ∇Wb(ra)⊗ (rb − ra) )−1 . (5.18) O uso desta técnica permite que o campo de velocidades seja corretamente avaliado [2]. 5.4 Movimento das Partículas Para um movimento mais ordenado das partículas pode ser utilizado o XSPH, que corrige a velocidade das partículas, é dado por: dra dt = va + ε ∑ b mb τab vbaWab (5.19) onde ρ = 1 2 (ρa + ρb) e ε é uma constante, onde seus valores variam de zero até um, sendo ε = 0.5 frequentemente usado. A velocidade é recalculada tendo em conta a velocidade dessa partícula e a velocidade média de todas as partículas que interagem com a partícula a. Assim são evitados interpenetração das partículas e as partículas �cam organizadas como pode-se notar na Figura(5.1), melhorando a estabilidade numérica das simulações. Devido ao suporte compacto do núcleo, somente as partículas vizinhas serão incluídas. Figura 5.1: Escoamento de lava, sem correção XSPH (a) e com (b). [18] 5. Correções Aplicadas ao SPH 39 5.5 Correção da Tensão Na implementação do SPHysics, existe a preocupação nos casos em que há o agru- pamento de partículas, que é denominada instabilidade de tensão. Um correção desta instabilidade, proposta por [13] e usada no código, é mostrada a seguir. Em [20] é feito um estudo sobre critérios de estabilidade para equações de SPH, e observou-se as seguintes condições de crescimento instável:∑ W ′′(r− r′).T > 0 (5.20) onde W ′′ é a segunda derivada do núcleo e a tensão T é negativa na compressão e positiva na tensão, como observado na �gura (5.2). Figura 5.2: Estudo da tensão SPH. [20] Na �gura podemos observar os regimes de estabilidade para o núcleo de spline cúbico. Se a derivada segunda for positiva (W ′′ > 0), eles são instáveis em tensão (T > 0) e estáveis na compressão. Se a derivada secundária é negativa, é instável na compressão e estável na tensão. A instabilidade de tensão resulta em um agrupamento de partículas de SPH. O agru- pamento é claro em materiais com uma equação de estado que pode dar origem a pressões negativas. O aglomerado das partículas de SPH não é físico porque será impedido em um sólido real pelas forças repulsivas entre os átomos. Em [13] é mostrado como a instabilidade no caso de �uidos pode ser removida usando a pressão arti�cial. Então basta adicionar a pressão arti�cial na equação do momento, como: dva dt = − ∑ b mb ( Pa ρ2 a + Pb ρ2 b + ∏ ab +Rfnab ) ∇aWab + g (5.21) O termo adicionado é chamado de termo de correção de tensão. Considerando um núcleo, a força repulsiva deve aumentar e a separação entre duas partículas diminuir. Para remover 5. Correções Aplicadas ao SPH 40 a instabilidade numérica, esta força repulsiva é escrita em termos do núcleo. Uma função adequada que aumenta as reduções de separação é fab = W (q) W (∇p) (5.22) onde ∇p é o espaçamento médio das partículas dividido pelo comprimento de suavização h: usando q = ∇p→ W (q) = W (∇p)→ fab = 1 usando q > ∇p→ W (q) < W (∇p)→ fab → 0 usando q < ∇p→ W (q) > W (∇q)→ fab >> 1 (área de interesse a ser corrigida). 5.6 Atualização Temporal Tendo a �nalidade de desenvolver as equações SPH no tempo, é aconselhável usar um esquema preciso de integração, pelo menos segunda ordem no tempo. Assim é utilizado o algoritmo Preditor-Corretor descrito por [10]. Considere a Equação do momento (4.1), Equação de continuidade (4.4), posição (5.19) e densidade de energia (4.7) da seguinte forma: dva dt = Fa (5.23) dρ dt = Da (5.24) dra dt = Va (5.25) dea dt = Ea (5.26) onde Va representa a contribuição de velocidade da partícula a e de partículas vizinhas (correção XSPH). O esquema preditor-corretor prevê a evolução no tempo como: vn+1/2 a = vna + ∇t 2 Fna (5.27) ρn+1/2 a = ρna + ∇t 2 Dn a (5.28) rn+1/2 a = rna + ∇t 2 Vn a (5.29) en+1/2 a = ena + ∇t 2 En a (5.30) calculando P n+1/2 a = f(ρ n+1/2 a ) de acordo com a equação (4.2). Esses valores são corrigidos usando forças no meio passo vn+1/2 a = vna + ∇t 2 Fn+1/2 a (5.31) ρn+1/2 a = ρna + ∇t 2 Dn+1/2 a (5.32) rn+1/2 a = rna + ∇t 2 Vn+1/2 a (5.33) 5. Correções Aplicadas ao SPH 41 en+1/2 a = ena + ∇t 2 En+1/2 a (5.34) Finalmente, os valores são calculados no �nal do intervalo de tempo seguinte vn+1 a = 2vn+1/2 a − vna (5.35) ρn+1 a = 2ρn+1/2 a − ρna (5.36) rn+1 a = 2rn+1/2 a − rna (5.37) en+1 a = 2en+1/2 a − ena . (5.38) Finalmente, a pressão é calculada a partir da densidade usando P n+1 a = f(ρn+1 a ). Segundo [10], esse desenvolvimento é usado para que o SPH conserve o momento linear e angular. Ná prática, eles usam o valor do ponto médio do tempo anterior em vez de calcular o valor no instante n, o que economiza tempo e cria apenas um pequeno erro. O esquema geral é de segunda ordem. 5.7 Simulador As simulações são efetuadas usando o SPHysics que é um código livre de hidrodinâ- mica de partículas suavizadas (SPH) inspirado na formulação de [11]. É uma colaboração conjunta entre vários pesquisadores da Universidade Johns Hopkins (EUA), da Univer- sidade de Vigo (Espanha), da Universidade de Manchester (U.K.) e da Universidade de Roma La Sapienza (Itália). O código SPHysics pode simular vários fenômenos, incluindo quebra de onda, barragens, objetos deslizantes, impacto de onda em uma estrutura etc. O modelo foi colocado em forma modular e uma variedade de recursos estão disponíveis para escolher diferentes opções de compilação. Figura 5.3: Árvore do diretório. O SPHysics é distribuído em um arquivo compactado (gz ou zip). A árvore de diretório é mostrada na Figura (5.3). Diferentes arquivos de saída são criados pela SPHysicsgen. Esses arquivos podem ser usados pelo executável SPHysics como arquivos de entrada ou por códigos MATLAB para visualizar os resultados (diferentes códigos MATLAB são fornecidos no subdiretório pós-processamento). Neste trabalho é usado um código em MATLAB que lê os arquivos de saída do SPHy- sics e permite gerar os grá�cos do campo de velocidades, pressão e posição das partículas. 5. Correções Aplicadas ao SPH 42 Assim podemos analisar os resultados de acordo com os parâmetros e as correções adota- das. Na �gura (5.4) são apresentados os passos de simulação, que serão aplicadas para solução da equação do momento, e as atualizações das propriedades do �uido, bem como a atualização das partículas no passo temporal. Figura 5.4: Esquema da simulação. Capítulo 6 Resultados A massa de �uido com uma con�guração retangular inicial H0 = 2m de altura e 1m de largura é deixada livre para �uir dentro de uma caixa de medidas 4m × 4m como na Figura (6.1). Figura 6.1: Con�guração inicial do tanque 6.1 Campo de velocidades A primeira análise é do campo de velocidades, onde é escolhido o núcleo Quíntico ou Wendland (3.24), e coe�ciente de viscosidade arti�cial α = 0.3, que é o sugerido pelas bibliogra�as, �ltro de densidade MLS (5.6), com correção a cada 30 passos, e velocidade do som c = 30m/s. Nas simulações foram utilizadas 5123 partículas. Nas �guras (6.2) e (6.3) são mostrados diferentes instantes da evolução na barragem. A barra de cores é comun a todos os instantâneos. As distâncias estão em metros e as velocidades em metros por segundo. Cada partícula é representada por uma cor corres- pondente à sua velocidade instantânea. Faz se uma comparação das velocidades com e sem o uso da viscosidade arti�cial. Percebe-se que o �uido escoa mais rapidamente sem o uso da viscosidade arti�cial, alcançando a parede sólida do tanque em 0.67s como mostrado na �gura (6.2). Já com a viscosidade arti�cial a parede sólida do tanque é alcançada em 0.77s. Na �gura (6.3) veri�ca-se um desprendimento das partículas ao se chocar com a parede sólida do tanque, o que não ocorre quando usada a viscosidade arti�cial como observado na mesma �gura. O mesmo desprendimento de partículas foi reportado em [7]. O com- portamento do campo de velocidades com o uso da viscosidade arti�cial está de acordo com as bibliogra�as pesquisadas, como [17]. 43 6. Resultados 44 Figura 6.2: Campo de velocidades do �uidos parte 1. Figura 6.3: Campo de velocidades dos �uidos parte 2. Por �m na �gura (6.3) em t=3s no caso sem uso da viscosidade arti�cial, há algumas partículas totalmente desprendidas, estas mantêm suas propriedades físicas inalteradas (pela cor da partícula), em relação a t=1.21s. Isto acontece porque estas partículas não 6. Resultados 45 tinham vizinhas em seu núcleo, logo suas propriedades físicas (nesse caso a velocidade) mudam a partir do momento que esta partícula possuir outras vizinhas. Isto con�rma a in�uência das partículas vizinhas no SPH, a�m de modelar propriedades físicas do �uido. 6.2 Campo de pressões O modelo tradicional do SPH não consegue recuperar o campo de pressões de forma suave, faz-se necessário o uso de correções para este �m, como �ltro de densidade. A seguir é apresentado a análise com respeito ao campo de pressões. Nas simulações foram utilizados 5123 partículas. A análise para o campo de pressões é feita com a �nalidade de avaliar o uso dos �ltros de densidade de Shepard e MLS. Novamente é escolhido o núcleo de Quíntico ou Wendland (3.24), e α = 0.3 para viscosidade arti�cial. Na �gura (6.4) são dispostos os comportamentos dos valores do campo de pressões para cada tempo, correspondentes ao escoamento simulado, para t=0.49s, 1.59s e 3.09s . É possível notar que o uso das correções Shepard e MLS conseguem prever o campo de pressões de forma suavizada. Através da �gura (6.5) nota-se a diferença das variações dos valores das pressões, referentes as simulações presentes na �gura (6.4). Os �ltros de densidade Shepard e MLS mantêm praticamente os mesmos valores, sendo o MLS com uma pequena vantagem, como se observa na Tabela (6.4). Figura 6.4: Pressão com uso da viscosidade arti�cial α = 0.3, com �ltros de densidade de Shepard e MLS. 6. Resultados 46 Figura 6.5: Grá�co dos valores da pressão com uso da viscosidade arti�cial α = 0.3 e com �ltros de densidade de Shepard e MLS. Tabela 6.1: Valores médios da pressão de acordo com o tempo, referentes as análises da �gura (6.5) Tempo Média das pressões Com viscosidade arti�cial e com �ltro de densidade de Shepard 0.49 6.1586× 103 1.59 1.9805× 103 3.09 2.3824× 103 Com viscosidade arti�cial e com �ltro de densidade MLS 0.49 6.1136× 103 1.59 2.0023× 103 3.09 2.3108× 103 6.3 Frente de onda: Validação experimental, tanque 1. Em [4] é desenvolvido um experimento que fornece dados da evolução da frente de onda com propriedades de �uidos diferentes. Dentre estes é utilizado a água, com a geometria descrita na �gura (6.6). As �guras (6.8) e (6.9) apresentam as comparações das simulações, usando o SPHysics com os resultados experimentais de [4], usa-se a seguinte equação para o cálculo das frentes de onda [7]: 6. Resultados 47 Figura 6.6: Geometria do tanque para validação experimental. [4] M xp = ( xexp − xSPH xexp ) (6.1) onde: M xp é a diferença percentual entre as abcissas das frentes de onda, obtidas experimen- talmente e pelo SPH, xexp é abcissa da frente de onda obtida experimentalmente, xSPH é abcissa da frente de onda obtida pelo SPH. Figura 6.7: Geometria 1 (a) e geometria 2 (b) para simulações. Na �gura (6.7) está disposta a con�guração inicial dos tanques para as simulações, com a �nalidade de comparar com os resultados experimentais. A primeira geometria a esquerda contém 3369 partículas e a segunda à direita 5854. Pelos resultados expostos na tabela (6.3) podemos constatar que a diferença entre as frentes de onda da simulação 6. Resultados 48 e com os resultados experimentais é bem pequena. Observa-se que com a evolução do passo temporal o erro percentual decai. Quando utilizado α = 0.3 foi possível alcançar melhores resultados. Tabela 6.2: Resultados das �guras (6.8) e (6.9), primeira geometria. Tempo (s) xexp(cm) xSPH Erro Percentual α = 0.2 0.1 19.0 17.7 6.84 0.2 29.0 29.7 -2.41 α = 0.3 0.1 19.0 17.5 7.89 0.2 29.0 29.1 -0.34 Figura 6.8: Comparação com modelo experimental, para α = 0.2 . Tabela 6.3: Resultados das �guras (6.10) e (6.11), segunda geometria. Tempo (s) xexp(cm) xSPH Erro Percentual α = 0.2 0.1 22.0 19.1 13.18 0.2 36.0 36.2 0.05 α = 0.3 0.1 22.0 18.5 15.90 0.2 36.0 34.9 3.00 6. Resultados 49 Figura 6.9: Comparação com modelo experimental, para α = 0.3 Figura 6.10: Comparação com modelo experimental, segunda geometria, para α = 0.2 6.4 Frente de onda: Validação experimental, tanque 2. Recentemente em [9] foi desenvolvido um experimento que fornece dados da evolução da interface de �uido; também neste caso o �uido é a água. A geometria e o dispositivo utilizado para o experimento, são descritos na �gura (6.12). 6. Resultados 50 Figura 6.11: Comparação com modelo experimental, segunda geometria, para α = 0.3 Figura 6.12: Dimensões do tanque 2 É utilizado a equação (6.1) para o cálculo da diferença entre a frente de onda modelada e do experimental. A tabela (6.4) mostra os resultados numéricos com uso da equação (6.1), concluindo que o modelo simulado SPH com uso da viscosidade arti�cial α = 0.1 obteve êxito. Foram utilizadas nessa simulação 36823 partículas. 6. Resultados 51 Figura 6.13: Comparação da interface do �uido simulada com experimento laboratorial, para os tempos 0.01s e 0.16s Figura 6.14: Comparação da interface do �uido simulada com experimento laboratorial, para os tempos 0.27s e 0.37s Figura 6.15: Comparação da interface do �uido simulada com experimento laboratorial, para o tempo 0.45s Nas �guras (6.13), (6.14) e (6.15) são expostas gra�camente as interfaces do �uido simuladas e comparadas com as correspondentes no experimento apresentado em [9]. 6. Resultados 52 Tabela 6.4: Resultados das �guras (6.13) e (6.14), segunda geometria. Tempo (s) xexp(cm) xSPH Erro Percentual α = 0.1 0.16 86.5 82.7 4.39 0.27 114.0 110.0 3.50 0.37 136.2 139.4 -2.34 Capítulo 7 Conclusão Os estudos realizados durante o projeto de mestrado forneceram subsídios para a mode- lagem do problema ruptura de barragem através do software SPHysics (open source), com a metodologia SPH. A formação nessa área do conhecimento se mostra muito importante, tendo em vista que o método estudado possui várias áreas de aplicações, possibilitando diversos estudos durante muitos anos de pesquisa. Nas apresentações feitas em eventos Locais (II Seminário de Métodos numéricos em Engenharia, 2017), Internacional (Colóquio Brasileiro de Matemática - IMPA 2017) e Nacional (Congresso Nacional de Matemática Aplicada e Computacional 2017) observou- se que o estudo do método SPH tem sido bem aceito e isto motiva a continuação dos estudos nessa temática. O SPHysics é um código aberto que se mostrou muito útil para a simulação do pro- blema da ruptura de barragem, além de possuir uma gama de opções que ainda podem ser exploradas. Foi possível analisar o uso da viscosidade arti�cial e observar que esta é uma ferramenta muito importante para o escoamento do �uido, conforme pode ser visto nas comparações feitas com o modelo experimental e também para a correção do campo de valocidades e pressão. O uso dos �ltros de Shepard e MLS também ajudaram a melhorar a simulaçao do campo de pressões. Além disso constatou-se, como já previsto nos estudos teóricos, que os �ltros de densi- dade de Shepard e MLS desempenham um importante papel no campo de pressões, sendo que este último possui desempenho teórico um pouco melhor devido sua capacidade de recuperar campos lineares. Vale ressaltar que os �ltros MLS e Shepard são usados pra resolver o mesmo tipo de di�culdade, mas em alguns casos o �ltro de Shepard pode ser considerado mais adequado pois o MLS é mais oneroso computacionalmente. A partir da experiência obtida será possível, em trabalhos futuros, explorar outros casos que o SPHysics possa nos ajudar. Ainda, para o problema de ruptura de barragem, pode-se modelar a ruptura de barragem no caso 3D, ou, ainda, no caso 2D considerar obstáculos no escoamento do �uido. Outra possibilidade de futuro vai na direção de pro- blemas com geometria mais exigentes, bem como casos multifásicos. Nesta direção grupos internacionais pesquisam ativamente, por exemplo, desenvolvendo o software DualPhysics que já inclui paralelização que aumenta signi�cativamente a performance nas simulações. 53 Referências [1] H.R. Alves. O rompimento de barragens no brasil e no mundo: Desastres mistos ou tecnolÓgicos? DOM, 2015. [2] T.S.L. Bonet; J. and Lok;. Variational and momentum preservation aspects of smo- othed particle hydrodynamics formulations. Computat.Methods Appl. Mech. Engine- ering, 180:97�115, 1999. [3] A.J. Crespo. Application of the Smoothed Particle Hydrodynamics model SPHysics to free-surface hydrodynamics. PhD thesis, Universidade de Vigo, Vigo-ESP, Junho 2008. [4] E.T. Cruchaga; M.A. Celentano; D.J. Tayfun. Collapse of a liquid column: numerical simulation and experimental validation. Comput. Mech., 39:453�476, 2007. [5] J.A. Cuminato; M. Meneguette. Discretização de Equações Diferenciais Parciais: técnicas de diferenças �nitas. SBM, 2013. [6] R.A. Fonseca. Algoritmos E�cientes em Métodos sem Malha. PhD thesis, UFMG, Minas Gerais, Março 2011. [7] F.D.C. Fraga. Estudo da fase gravitacional-inercial do espalhamento de óleo em mar calmo empregando o método lagrangiano de partículas Smoothed Particle Hydrody- namics. PhD thesis, UFES, Espírito Santo, Dezembro 2014. [8] R.G. Liu; M.B. Liu. Smoothed Particle Hydrodynamics, a meshfree particle method. World, Scienti�c Publishing, 2003. [9] A. Lobovsky; L. Botia-Vera; E. Castellana; F. Mas-Soler; J. Souto-Iglesias. Experi- mental investigation of dynamic pressure loads during dam break. Journal of Fluids and Strutures, 48:407�434, 2014. [10] J.J. Monaghan. On the problem of penetration in particle methods. Journal Com- putational Physics, 82:1�15, 1989. [11] J.J. Monaghan. Smoothed particle hydrodynamics. Annual Rev. Astron. Appl., 30:543�574, 1992. [12] J.J. Monaghan. Simulating free surface �ows with sph. Journal of Computational Physics, 110:399�406, 1994. [13] J.J. Monaghan. Sph without tensile instability. Journal Computational Physics, 159:290�311, 2000. [14] J.J. Monaghan. Conduction modelling using smoothed particle hydrodynamics. Re- ports on Progress in Physics, 68:1703�1759, 2005. 55 REFERÊNCIAS 56 [15] J.J. Monaghan and Kos A. Solitary waves on a cretan beach. Journal Waterway, Portm Coastal and Ocean Engineering., 125:145�154, 1999. [16] J.P. Morris; P.J. Fox; Y. Zhu. Modeling low reynolds number incompressible �ows using sph. Journal of Computational Physics, 136:214�226, 1997. [17] E.A.P. Narino. Modelagem numérica de micro�uídica através do Método Lagrangeano sem malha Smothed Particle Hydrodynamics. PhD thesis, Unicamp, Campinas, 2017. [18] A. Paiva. Uma abordagem lagrangeana para simulação de escoamentos de �uidos viscoplásticos e multifásicos. PhD thesis, PUC-RJ, Rio de Janeiro, 2007. [19] A. Paiva; F. Petronetto; G. Tavares; T. Lewiner. Simulação de Fluidos sem Malha: Uma introdução ao método SPH. IMPA, 2009. [20] W.J. Swegle. Smoothed particle hydrodynamics stability analysis. Journal of Com- putational Physics, 116:123�134, 1995. [21] J.R.G. Vasco. Desenvolvimento de software utiliando a técnica SPH na geração de ondas de submersão. PhD thesis, Unesp, 2014. [22] B.D. Violeau; D. Rogers. Smoothed particle hydrodynamics (sph) for free-surface �ows: past, present and future. Journal of Hydraulic Research, 54:1�26, 2016. [23] R. Xu; P. Stansby; D. Laurence. Accuracy and stability in incompressible sph (isph) based on the projection method and a new approach. Journal of Computational Physics, 228:6703�6725, 2009.