Joel Roberto Guimarães Vasco Desenvolvimento de software utilizando a técnica SPH (Smoothed Particle Hydrodynamics) na geração de ondas de submersão Ilha Solteira 2014 Joel Roberto Guimarães Vasco Desenvolvimento de software utilizando a técnica SPH (Smoothed Particle Hydrodynamics) na geração de ondas de submersão Tese de Doutorado apresentada à Faculdade de Engenharia de Ilha Solteira – Unesp – como parte dos requisitos para obtenção do título de Doutor em Engenharia Elétrica. Especialidade: Automação Orientador: Carlos Roberto Minussi Coorientador: Geraldo de Freitas Maciel Ilha Solteira 2014 Vasco Desenvolvimento de soIlha Solteira2014 124 Sim Tese (doutoEngenharia AutomaçãoSim . . FICHA CATALOGRÁFICA Desenvolvido pelo Serviço Técnico de Biblioteca e Documentação Vasco, Joel Roberto Guimarães Vasco. Desenvolvimento de software utilizando a técnica sph (smoothed particle hydrodynamics) na geração de ondas de submersão / Joel Roberto Guimarães Vasco Vasco. -- Ilha Solteira: [s.n.], 2014 124 f. : il. Tese (doutorado) - Universidade Estadual Paulista. Faculdade de Engenharia de Ilha Solteira. Área de conhecimento: Automação, 2014 Orientador: Carlos Roberto Minussi Co-orientador: Geraldo De Freitas Maciel Inclui bibliografia 1. Geração de ondas. 2. Métodos numéricos. 3. Ondas de impacto. 4. Métodos particulados. 5. SPH. 1,2,3,4,5 = 690 V331d Para meus pais AGRADECIMENTOS Mais uma etapa de um longo caminho encontra seu fim. Neste difícil percurso pessoal, somos agraciados com o convívio de várias pessoas, que também fazem suas jornadas e em algum momento há uma interseção de caminhos. Nestas interseções, auxílios são naturalmente ofertados, mesmo na forma de uma conversa informal, e são capazes de garantir a energia e disposição necessárias para continuarmos a jornada. É neste espaço oportuno que essas pes- soas (ou entidades) são lembradas. Listarei todos aqueles que de alguma forma contribuíram para a conclusão com sucesso do trabalho em questão, mas não por ordem de importância. De antemão, peço desculpas a alguém que porventura venha a esquecer. Isto posto, agradeço: • À Agência de fomento FAPESP pela concessão da bolsa, garantindo minha dedicação exclusiva à pesquisa. Agradeço também ao Programa de Pós-Graduação em Engenharia Elétrica (PPGEE), por entender e abraçar meu plano de pesquisa. • À CAPES/FCT, que através do projeto “Amigos de Boussinesq”, do qual faço parte há muito tempo, possibilitou-me um período sanduíche de um ano em Lisboa-PT. Neste con- texto, não posso deixar de citar os pesquisadores do LNEC que tão bem me receberam: Juana Fortes, Eric Didier, Rui Capitão, Ruth Lemos, João Alfredo e Teresa Reis. Tam- bém agradeço ao meu tutor, Luís Gil, da Universidade Nova de Lisboa, que em muito contribuiu para garantir uma boa estada em terras além-mar. • Aos colegas de trabalho do LNEC: Jorge Gadelho, Maria de Fátima, Gonçalo Jesus, Ru- ben Silva, Diogo Neves, Gil Souza, Euclides Rodrigues, Ricardo Martins, Sara Rodrigues, Edgar Matias, Ana Mendonça, Laurinda Almada, Mariana Rocha, Joanne Laranjeiro, Má- rio Moreno, Gonçalo Faria, Javier Ramalleira e Liliana Pinheiro, pela convivência. • Ao Laboratório do Tanque de Provas Numérico (TPN) da Poli-USP, na figura dos profes- sores Cheng Liang Yee e Kazuo Nishimoto, que sempre foram extremamente receptivos em minhas visitas. Uma pena que as visitas ficaram, com o passar dos anos, menos fre- quentes do que eu gostaria. Aos colegas do TPN: Márcio Tsukamoto, Fábio Motezuki e Pedro Cardozo. • À banca examinadora, por aceitar contribuir, com entusiasmo e critério científico, na ava- liação deste trabalho. • Aos docentes de Ilha Solteira, pelo convívio de longa data. Agradeço aos meus orien- tadores, na figura do Prof. Carlos Minussi, pelo meu acolhimento do projeto dentro do programa de pós-graduação, além de ser uma pessoa com atitude extremamente otimista, que propõe soluções (concretas) aos problemas vivenciados, além das oportunidades vis- lumbradas e concretizadas. Agradeço ao Prof. Milton Dall’Aglio, pela capacidade em resolver problemas em diversos assuntos e bons conselhos. Apenas lamento a falta de tempo para poder aprender como ele tudo o que julgo suficiente. Agradeço também ao Prof. Dib Gebara, pelas conversas e conselhos. • Aos colegas de trabalho de Goiânia, pelo acolhimento e carinho com que fui recebido. • Aos colegas pessoais, quer em São Paulo, quer em Ilha Solteira: Humberto Ruggeri, Evandro Cunha, Luciana Kataoka, Adriana Vieira, Fabiana Oliveira, José Carlos, San- dra e Guilherme Fiorot. • À minha família, principalmente pelo apoio incondicional. Agradeço também aos meus amigos de minha terra natal, Orlando Cavalari e Fabrício de Paula. RESUMO O escorregamento de massas sólidas em lagos de águas tranquilas geram ondas devido à transferência de energia da massa deslizante para o corpo d’água. Essas ondas, chamadas de ondas de submersão, em particular as ondas solitárias, são conhecidas e estudadas há vários anos principalmente pela capacidade energética deste tipo de onda e seu potencial danoso. Entretanto, aproximações analíticas do fenômeno em tela esbarram em diversas dificulda- des. Em outras palavras, a fase de geração de ondas, notadamente das ondas de impacto, é sabidamente complexa, não apenas pela dificuldade no estabelecimento da superfície livre (região do splash), mas também pelo relacionamento da dinâmica do material impactante com a altura da onda resultante. Por esse motivo, técnicas numéricas vêm sendo emprega- das nestes fenômenos. Entre as aproximações numéricas utilizadas mais recentemente, os métodos sem malha vêm ganhando notoriedade, principalmente frente às suas vantagens. Livres da discretização da malha, problemas de superfície livre, mesmo com alta variabi- lidade, são factíveis. Neste contexto, lança-se mão de um código baseado no modelo SPH (Smoothed Particle Hydrodymanics), desenvolvido no âmbito desta Tese de Doutoramento e validado utilizando casos clássicos da literatura em escoamentos de fluido ideal (ruptura de barragens, geração e quebra de onda e esvaziamento de um reservatório) e real, de reo- logia newtoniana (Poiseuille com superfície livre, Poiseuille plano, escoamento de Couette e ruptura de barragens). De posse do código numérico, estima-se a altura da onda solitária gerada a partir do impacto de uma massa deslizante em meio líquido. A massa deslizante é representada por duas formas distintas: um bloco indeformável, modelado como um corpo rígido e um fluido ideal. No caso do bloco indeformável, embora o perfil de onda seja bem representado, a altura da onda é superestimada, fato também relatado pela literatura no assunto. No caso da massa deslizante modelada por um fluido ideal, há concordância numérico-experimental tanto na forma quanto na amplitude da onda gerada. Palavras-chave: Onda solitária. Impacto. SPH. ABSTRACT Landslides that impact into water (such as lakes) can generate waves due to the energy transfer from the solid to the water. These submerged waves, or more specifically solitons, are well known and have been studied for many years specially because of their massive energy and destructive potential. However, analitical approches of this problem are difficult. In another words, the wave’s generation phase is complex, primary due to the determination of the free surface (splash) and correlation between the energy of the landslide and the genereted wave. So, numerical methods are an alternative and interesting option. Among the recently used numerical methods, the mesh free has become notorious because of its advantages. Without the mesh restrain, free surface problems are easier to handle. In this context it is proposed a program based on SPH (Smoothed Particle Hydrodymanics) model, developed in this Thesis and tested in different scenarios considering the flow of ideal (dam break, wave generation and wave breaking and reservoir) and newtonian fluids (free surface Poiseuille flow, plane Poiseuille flow, Couette flow and dam break). So, the numerical code is used in order to evaluate the solitary wave height that is generated from the impact of a mass, which represents a landslide in a fluid. The mass is represented by two forms: an undeformable solid block and an ideal fluid. Both forms are tested, and the solitary wave is reproduced, and compared to experiments. Differences are noted among wave heights using a solid block, as reported by the literature. The comparison with the results from the mass represented as an ideal fluid showed a good approach for both the shape and the wave amplitude. Keywords: Soliton. Impact. SPH. LISTA DE ILUSTRAÇÕES Figura 1 – Exemplo de problema complexo: uma comparação entre um experimento (à esquerda) e o modelo numérico SPH desenvolvido pelo autor (à direita), em instantes aproximadamente iguais. 34 Figura 2 – À esquerda, o raio de ação da função W e à direita os valores da grandeza interpolados, representados pelas barras verticais cuja altura é proporcional ao inverso da distância ri j. 42 Figura 3 – Subdivisão do domínio computacional para diminuir o número de partículas de busca na criação da lista de vizinhança. À esquerda, vê-se todo o domínio computacional, enquanto à direita vê-se o detalhe da região de interesse. 43 Figura 4 – Função de suavização ou núcleo W . 44 Figura 5 – Partículas fantasma, com condição �ugh ·�t =−�ui ·�t, ∀i. 51 Figura 6 – Croqui esquemático com as distâncias horizontais (ψ) e verticais (ξ ) entre as partículas da fronteira e fluida. 54 Figura 7 – À esquerda, um croqui esquemático de uma ruptura de barragens, enquanto que à direita, observa-se o líquido retido, em ensaio experimental. 59 Figura 8 – Comparação numérico-experimental, para uma relação H0/L0 = 1 com (a) H0 = 2,86 cm e (b) H0 = 5,72 cm. 61 Figura 9 – Comparação numérico-experimental, para uma relação H0/L0 = 2 com (a) H0 = 5,72 cm e (b) H0 = 2,86 cm. 61 Figura 10 – Disposição geral para geração de ondas: (a) croqui esquemático e (b) domí- nio discretizado. 63 Figura 11 – Comparação teórico-numérica da frequência da onda gerada pelo batedor. 64 Figura 12 – Comparação teórico-numérica da altura H/S da onda gerada pelo batedor. 66 Figura 13 – Posição da superfície livre, para o tempo t = 3,91 s, com w = 5,5 Hz, S = 0,25 m, declividade da praia de 10% e a profundidade normal de 0,5 m 67 Figura 14 – À esquerda, são apresentadas as características geométricas do reservatório e orifício simulados. À direita, percebe-se a disposição inicial das partículas para um nível de aproximadamente 1,2 m no reservatório. 68 Figura 15 – Comparação teórico-numérica da variação do nível do reservatório. 69 Figura 16 – Evolução temporal no problema de esvaziamento de um reservatório. As diferentes figuras (da esquerda para a direita, de cima para baixo), corres- pondem aos tempos 0,53, 1,38, 2,63 e 3,61 s de simulação. 70 Figura 17 – Função de suavização, derivada primeira e derivada segunda de W , consi- derada como sendo do tipo spline cúbica 2D normalizada neste exemplo. Como dito no texto, percebe-se uma mudança de sinal na derivada segunda com a distância adimensional s. (Nota: as derivadas são calculadas por mé- todos numéricos de primeira ordem). 73 Figura 18 – Perfil de velocidade para o Poiseuille plano permanente (à esquerda), para um Re ≈ 0,01 e resultados de Cueille (2005), para o Poiseuille plano, com Re = 0,01 e sem suavização do campo de velocidades (à direita). 75 Figura 19 – Comparação entre o perfil de velocidades esperado (linha contínua, teórica) com o numérico (pontos). Mesmo com a presença da viscosidade, as partí- culas se mantém organizadas, fato esse ilustrado pelos pontos na figura, que representam uma linha de partículas, que permanecem com a mesma altura do começo da simulação. 77 Figura 20 – À esquerda, a variação espacial da massa específica, no tempo t = 1,5 s, mostra que a condição de incompressibilidade é reproduzida no escoamento. À direita, a posição das partículas no mesmo instante. 78 Figura 21 – À esquerda, disposição das partículas nos instantes iniciais e, à direita, o escoamento no instante t = 0,29 s. O domínio é retratado apenas em um trecho de 1 m, sendo a fronteira periódica capaz de reproduzir, na simulação, o regime permanente. A coloração indica a intensidade do vetor velocidade (azul escuro, u = 0, vermelho, u = umax). 79 Figura 22 – À esquerda, comparação entre o perfil de velocidade analítico e numérico, com variação na resolução espacial. À direita, influência do tipo de veloci- dade atribuída às partículas fantasma, para condição de contorno de aderên- cia no fundo do canal, com d/Δx = 35. 79 Figura 23 – À esquerda, disposição das partículas nos instantes iniciais e, à direita, o escoamento no instante t = 0,20 s. O domínio é retratado apenas em um trecho de 1 m e a coloração indica a intensidade do vetor velocidade (azul escuro, u = 0, vermelho, u = umax). 80 Figura 24 – Resultados analíticos e numéricos do perfil de velocidade do Poiseuille plano no centro do domínio computacional , nos tempos t = 0,01 s (+), 0,02 s (×), 0,05 s (∗), 0,1 s (�) e ∞ (�). Os símbolos indicam os resultados numéri- cos, enquanto que a linha contínua representa os resultados analíticos, nos respectivos tempos. 81 Figura 25 – À esquerda, disposição das partículas nos instantes iniciais e, à direita, o escoamento no instante t = 0,10 s. O domínio é retratado apenas em um trecho de 1 m, sendo a fronteira periódica capaz de reproduzir, na simulação, o regime permanente. A coloração indica a intensidade do vetor velocidade (azul escuro, u = 0, vermelho, u = umax). 82 Figura 26 – Resultados analíticos e numéricos do perfil de velocidade do Couette no centro do domínio computacional, nos tempos t = 0,01 s (+), 0,02 s (×), 0,03 s (∗), 0,04 s (�) e ∞ (�). Os símbolos indicam os resultados numéri- cos, enquanto que a linha contínua representa os resultados analíticos, nos respectivos tempos. 83 Figura 27 – Comparação entre os resultados experimentais de Minussi (2007), teóricos de Ritter e numéricos para o problema de ruptura de barragem retendo ma- terial newtoniano, com altura de 0,07 m (a) e 0,10 m (b). 84 Figura 28 – Comparação entre os resultados experimentais de Leite (2009) e numéricos para o problema de ruptura de barragem retendo material newtoniano, com altura de 0,08 m (a), 0,10 m (b) e 0,12 m (c). 85 Figura 29 – Tipos de leis reológicas. 87 Figura 30 – Variação da viscosidade aparente com a profundidade, em um problema tipo Poiseuille com superfície livre, com μmax = 1000 Pa s. 88 Figura 31 – À esquerda, pode-se ver o canal de ondas onde os ensaios foram efetuados e à direita, um detalhe da rampa de lançamento de material deslizante com comporta aberta. 93 Figura 32 – Procedimento para determinação da velocidade da frente do deslizamento pelo tratamento de imagens sucessivas. 94 Figura 33 – Comparação numérico e experimental para a onda gerada pelo impacto de um bloco único 97 Figura 34 – Impacto do bloco e formação da onda, sendo as características da simulação as mesmas que as utilizadas em Monaghan, Kos e Issa (2003) 98 Figura 35 – À esquerda, bloco utilizado por Maciel e Nascimento (2002) e à direita, configuração da primeira aproximação. 99 Figura 36 – À esquerda, a comparação numérica para o cálculo de pressão em um tanque em repouso. À direita, a resultante do balanço de quantidade de movimento numérica gerada apenas pelo gradiente de pressão. O valor teórico de tal balanço deve ser igual a 9,81. 99 Figura 37 – Comparação numérico experimental da onda gerada pelo impacto de um bloco em águas tranquilas, em função do número de Froude no impacto. 100 Figura 38 – Geração numérica da onda solitária com o código SPH, utilizando as condi- ções de contorno dadas pelas forças radiais (MONAGHAN; KAJTAR, 2009). 101 Figura 39 – (a) Posição da superfície livre, medida a aproximadamente 1,7 m do ponto de impacto e (b) a dinâmica do centro de gravidade do material deslizante, para d = 0,15 m (8022 partículas). 103 Figura 40 – (a) Posição da superfície livre, medida a aproximadamente 1,7 m do ponto de impacto e (b) a dinâmica do centro de gravidade do material deslizante, para uma lâmina d’água de 0,175 m (6878 partículas). 103 Figura 41 – (a) Posição da superfície livre, medida a aproximadamente 1,7 m do ponto de impacto e (b) a dinâmica do centro de gravidade do material deslizante, para uma lâmina normal de 0,20 m (6035 partículas). 104 Figura 42 – Imagens para vários tempos (de cima para baixo): 16/30, 24/30, 33/30, 42/30 e 53/30 s. À esquerda, resultados experimentais e à direita as simu- lações SPH. 105 Figura 43 – À esquerda, resultado exato para interpolação proposta e à direita os resul- tados do interpolador MLS. 122 Figura 44 – À esquerda, resultado do SPH para interpolação da função f (x,y) = x+ y e à direita os erros envolvidos. 122 Figura 45 – À esquerda, resultado exato para interpolação proposta e à direita os resul- tados do interpolador MLS, agora com as partículas desordenadas. 123 Figura 46 – À esquerda, resultado do SPH para interpolação da função f (x,y) = x+ y e à direita os erros envolvidos, com as partículas desordenadas. 124 LISTA DE SÍMBOLOS E ABREVIAÇÕES α , β parâmetros da viscosidade artificial no método SPH βm fator de correção que garante que a força exercida pela parede nas partículas fluidas (equação 66, p. 54) seja independente do número de partículas da fronteira Δh incremento dos métodos de integração de equações diferenciais de primeira ordem (s) Δp espaçamento entre partículas na fronteira (m) Δt passo de tempo numérico (s) Δx espaçamento médio entre partículas (m) Δ incremento no raio de vizinhança para cálculo da tabela Verlet (m) δ (�r−�r j) função impulso ou delta de Dirac κ fator de proporcionalidade entre rs e Δx λ parâmetro do método MPS (equação 9, p. 38) 〈·〉 indica aproximação da função original |�D | taxa de deformação, calculado usando |�D |= √ �D : �D (s−1) �D tensor taxa de deformação (s−1) �T tensor de tensões (Pa) K índice de consistência (Pa s) N índice do escoamento Tc tensão crítica (Pa) �F segundo membro da equação da conservação de quantidade de movimento, escrita na forma SPH (m/s2) D segundo membro da equação da conservação de massa, escrita na forma SPH (kg/m3 s) H segundo membro da equação 34 (p. 44), quando o comprimento de suavização é consi- derado variável na simulação (m/s) μ viscosidade dinâmica do fluido (Pa s) μapp viscosidade aparente (Pa s) ν viscosidade cinemática do fluido (m2/s) φ campo qualquer, seja escalar, vetorial ou tensorial Πi j viscosidade artificial ψ distância horizontal entre as partícula de fronteira e fluida (ver fig. 6, p. 54, m) ρ massa específica do fluido (kg/m3) ρ0 massa específica inicial (kg/m3) ρs massa específica do material sólido (kg/m3) θ inclinação, seja de um canal ou rampa (rad) ε constante XSPH que varia entre 0 e 1 (valor usual: 0,5, equação 40, p. 45) ϖ função peso no método MPS (equação 4, p. 37) �∇ vetor nabla, de coordenadas ( ∂ ∂x , ∂ ∂y ) �Ω velocidade angular do material único (rad/s) �τ torque total atuando no corpo único (N m) �C centro de massa do material deslizante único (m, m) �F força total atuante no material único (N) �f força por unidade de massa atuando na partícula (m/s2) �fp função com as mesmas características do núcleo de suavização (MONAGHAN; KOS; ISSA, 2003) �j vetor unindo as partículas i e j �U velocidade do centro de massa do material único (m/s) �σ tensor de Cauchy (�σ =−p�I + �T ) (Pa) �g representa a força de corpo gravitacional (m/s2) �n versor normal externo à parede �r posição da partícula (x, y) (m, m) �t versor normal tangencial à parede �u vetor velocidade do escoamento (m/s) ξ distância vertical entre as partícula de fronteira e fluida (ver fig. 6, p. 54, m) ζi variáveis do interpolador MLS Ao área do orifício (vide fig. 14, p. 68, m2) Ar área do reservatório (vide fig. 14, p. 68, m2) Ai j uma matriz, geralmente relacionada a problemas lineares tipo Ax = b Bi j matriz de renormalização c celeridade do som (m/s) Cd coeficiente de descarga do orifício D dimensões do espaço de simulação (D = 2) d profundidade normal (m) e1, e2 parâmetros para cálculo da correção da instabilidade de tensão (equação 67, p. 55), que variam conforme pi e p j, respectivamente. Fx componente horizontal da força de corpo, no problema do Poiseuille plano Frimp número de Froude no impacto do bloco com o meio líquido G índice (sobrescrito) que representa o núcleo gaussiano gh índice (subscrito) que indica partícula fantasma H altura da onda (m) h comprimento de suavização (m) H0 altura inicial do fluido, utilizada em problemas do tipo ruptura de barragens ou no esva- ziamento de um reservatório (m) h0 comprimento de suavização inicial na simulação (m) I matriz identidade i, l índice (subscrito) para denominar partículas fluidas IS índice (sobrescrito) para indicar a aplicação da interpolação de Shepard J momento de inércia em torno do centro de massa, para o material único (kg m2) j índice (subscrito) para determinar partículas fluidas vizinhas da partícula i K uma constante k índice (subscrito) que indica partícula de parede k0 número de onda (m−1) L comprimento de onda (m) L0 largura inicial do material retido, em problemas do tipo ruptura de barragens (m) M massa do material deslizante único (kg) m massa (kg) Ma número de Mach MLS índice (sobrescrito) para indicar a aplicação da interpolação MLS n contador incremental N0 indicativo de massa específica inicial P expressão que determina a componente horizontal da força por unidade de massa exer- cida pela fronteira p pressão (Pa) p1, p2 constantes da equação 64 (p. 53, Monaghan, 1994) Q expressão que determina a componente vertical da força por unidade de massa exercida pela fronteira (m/s2) rp raio da partícula (m) rs domínio de suporte (m) ri j distância entre as partículas i e j, ou seja: ri j = ‖�ri −�r j‖ S deslocamento (stroke) do gerador de ondas (m) s mudança de variáveis para cálculo do núcleo de suavização: s = r/h SC índice (sobrescrito) que representa o núcleo dado por uma spline cúbica T tempo adimensional, dado por T = t √ H0 g/L2 0 t tempo (s) ux(y, t) Perfil de velocidades teórico, nos problemas do Poiseuille e Couette (m/s) V máxima velocidade encontrada na simulação (m/s) Vimp velocidade de impacto do bloco com o meio líquido, quando o centro de massa do bloco intercepta a lâmina d’água normal (m/s) vsig velocidade de referência (m/s) para cálculo da viscosidade artificial (KAJTAR; MO- NAGHAN, 2008) W representa a função peso, ou núcleo de suavização w frequência da onda (s−1) x direção ou coordenada horizontal (m) y direção ou coordenada vertical (m) Z alcance adimensional da frente, em problemas de ruptura de barragens z alcance da frente de onda, em problemas do tipo ruptura de barragens (m) MPS Moving Particle Semi-implicit Method SPH Smoothed Particle Hydrodymanics ALE Arbitrary Lagrange Euler EFGM Element-free Garlekin Method ISPH Incompressible Smoothed Particle Hydrodynamics LES Large Eddy Simulation MDF Método das Diferenças Finitas MEC Método dos Elementos de Contorno MEF Método dos Elementos Finitos MFI Método de Fronteira Imersa, ou Immersed Boundary Method MLS Moving Least Squares MSM Métodos Sem Malha MVF Método dos Volumes Finitos PIV Particle Image Velocimetry VOF Volume of fluid SUMÁRIO 1 INTRODUÇÃO 27 1.1 OBJETIVO 29 1.2 ORGANIZAÇÃO DO TRABALHO 29 2 REVISÃO BIBLIOGRÁFICA 31 2.1 FORMULAÇÃO EULERIANA × LAGRANGIANA 31 2.2 MÉTODOS NUMÉRICOS TRADICIONAIS 31 2.3 DIFERENÇAS ENTRE MÉTODOS TRADICIONAIS E MÉTODOS SEM MA- LHA (MSM) 33 2.4 PASSOS GERAIS PARA RESOLVER UM PROBLEMA USANDO MSM 35 2.5 PRINCÍPIOS DOS DIFERENTES MÉTODOS NUMÉRICOS 36 3 O MÉTODO MPS 37 4 O MÉTODO SPH 39 4.1 FUNDAMENTOS DO MÉTODO SPH 39 4.2 A FUNÇÃO NÚCLEO DE SUAVIZAÇÃO W 41 4.3 COMPRIMENTO DE SUAVIZAÇÃO 43 4.4 EQUAÇÃO DA QUANTIDADE DE MOVIMENTO 44 4.5 DESLOCAMENTO DE PARTÍCULAS 45 4.6 EQUAÇÃO DA CONTINUIDADE 45 4.7 CORREÇÕES APLICADAS AO SPH 46 4.7.1 Reinicialização da massa específica 47 4.7.1.1 Interpolador de Shepard 47 4.7.1.2 Interpolador MLS 47 4.7.2 Renormalização 48 4.8 LEIS DE ESTADO PARA PRESSÃO 49 4.9 VISCOSIDADE ARTIFICIAL 50 4.10 CONDIÇÕES DE CONTORNO 50 4.10.1 Partículas fantasma 51 4.10.2 Fronteiras reativas 53 4.11 INSTABILIDADES NO SPH 54 4.12 EVOLUÇÃO TEMPORAL 55 4.12.1 Preditor-Corretor 56 4.12.2 Método Verlet 56 4.12.3 Runge-Kutta 57 4.13 COMPARAÇÃO ENTRE MÉTODOS NUMÉRICOS: SPH e MPS 58 5 ESCOAMENTO DE FLUIDOS IDEAIS 59 5.1 RUPTURA DE BARRAGENS 59 5.1.1 Descrição do modelo numérico 60 5.1.2 Resultados 60 5.2 GERAÇÃO, PROPAGAÇÃO E QUEBRA DE ONDAS 62 5.2.1 Fase de geração 62 5.2.1.1 Teoria de pequena amplitude 63 5.2.1.2 Modelo numérico 63 5.2.1.3 Resultados 64 5.2.1.4 Teoria do gerador de onda 65 5.2.2 Propagação de ondas 65 5.2.3 Quebra de ondas 66 5.3 ESVAZIAMENTO DE UM RESERVATÓRIO 67 6 ESCOAMENTO DE FLUIDOS REAIS 71 6.1 ESCOAMENTO DE FLUIDOS NEWTONIANOS 71 6.1.1 Incorporação do tensor de tensões 72 6.1.2 Viscosidade numérica como viscosidade real 72 6.1.3 Demais aproximações 74 6.1.4 Inconvenientes da formulação adotada 75 6.2 ESTUDOS DE CASO 76 6.2.1 Fronteiras periódicas 76 6.2.2 Poiseuille com superfície livre 78 6.2.3 Poiseuille plano 80 6.2.4 Couette 81 6.2.5 Ruptura de barragem 82 6.3 ESCOAMENTO DE FLUIDOS NÃO NEWTONIANOS 84 6.3.1 Reologia 86 6.3.2 Tensor de tensões 86 7 ONDAS DE IMPACTO 89 7.1 REVISÃO BIBLIOGRÁFICA 90 7.2 METODOLOGIA 92 7.2.1 Procedimento experimental: material fragmentado 93 7.2.2 Procedimento experimental: material único 94 7.2.3 Procedimento numérico: dinâmica de um corpo rígido 95 7.2.4 Procedimento numérico: material fragmentado 96 7.2.5 Resultados 96 27 1 INTRODUÇÃO As ondas, sejam mecânicas ou eletromagnéticas, representam ao mesmo tempo um fenômeno complexo e de fundamental importância para toda humanidade. Tal assertiva pode ser verificada observando a gama de situações que tem ligação direta com o conceito de ondas, por exemplo: a propagação da luz natural do sol, a capacidade de distinguir sons diversos ou a quebra de ondas do mar na praia. As situações evocadas dão apenas uma ideia dos diferentes ramos em que as ondas são o fenômeno característico e principal objeto de estudo. Dada a importância e a frequência com que as ondas são observadas na natureza, os pes- quisadores vêm dispendendo esforços no sentido de compreender e evitar os problemas relacio- nados a esses fenômenos. Pode-se citar, com intuito ilustrativo, os prejuízos e as consequências que podem ser ocasionadas pelas ondas, em situações não controladas: a) prejuízos financeiros, através de uma falha no sistema de distribuição de energia elétrica; b) danos ambientais, causados por uma ruptura de uma barragem retendo rejeitos de indús- trias de papel; c) prejuízos materiais, através da ruína de estruturas portuárias devido à excessiva ondulação causada por mau tempo e; d) perda de vidas, como no caso de uma tsunami atingindo uma comunidade junto à linha de costa. Sendo assim, desenvolvem-se meios de evitar ou ao menos aliviar os efeitos danosos apontados, e que passam, invariavelmente, pela análise do padrão de onda. Não é preciso lem- brar, todavia, que as análises do padrão de onda só são possíveis graças a etapa de monitora- mento, o que muitas vezes ocasiona dificuldades por causa da quantidade excessiva de dados a serem tratados, tornando a intervenção uma tarefa mais lenta (ALLEN; APOSTOLV; KREISS, 2005). No caso do sistema de distribuição de energia elétrica, há uma necessidade em se obter, de forma rápida, o tipo de anomalia presente na onda de tensão, tais como transitórios, variação de tensão, distorções de forma, flutuações ou oscilações (MALANGE, 2010). A rapidez é ne- cessária para que ações de mitigação possam ser empreendidas, afim de evitar, em um caso mais extremo, um colapso no sistema (vulgo “apagão”). As ações de intervenção são balizadas por um tratamento prévio dos parâmetros da onda (frequência, amplitude, simetria e formato), obti- dos com a utilização de métodos tais como: Tranformadas de Fourier, Wavelets, Redes Neurais, Lógica Fuzzy ou uma combinação entre eles (MALANGE, 2010). 28 Capítulo 1. INTRODUÇÃO De forma semelhante, sistemas de previsão e alerta de navegação de embarcações ba- seiam-se no monitoramento dos ventos incidentes e principalmente ondas, e no tratamento dos dados. Tendo em vista o caráter estratégico das hidrovias no desenvolvimento do país, principal- mente no tocante ao transporte de mercadorias, projetos como o FINEP Ondisa, desenvolvido na Unesp - Ilha Solteira, visam integrar um sistema de monitoramento e previsão de altura de ondas para identificar rotas mais seguras para embarcações na hidrovia Tietê-Paraná. Não é pre- ciso ressaltar que o sucesso do sistema de alerta passa pelos modelos empregados no tratamento dos dados. Sendo assim, sabe-se que diversas técnicas numéricas e de programação foram desen- volvidas para lidar com o tratamento dos sinais, ou de maneira mais específica, da onda inci- dente. Por exemplo, Allen, Apostolv e Kreiss (2005) ressaltam os benefícios na substituição de um tratamento prévio dos dados de entrada através de programas específicos, em detrimento a dispositivos físicos, resultando em uma redução do tempo para a tomada das ações de mitiga- ção. Nesse sentido, alguns esforços também são dirigidos não só para o tratamento do sinal em si, mas também da sua geração. Uma onda pode ser gerada das mais variadas formas, sendo caracterizada por uma per- turbação ou pulso que tem capacidade de transportar energia. Por exemplo, podem-se citar as ondas geradas pelo deslocamento de uma embarcação em meio líquido (ondas de Kelvin, ou ship waves), ondas geradas pelo impacto da estrutura do casco de embarcações na superfície da água (slamming), ondas geradas por impacto de massas sólidas em meio líquido (deslizamento de terra ou rocha em lagos de barragem), entre outras. Modelos numéricos são desenvolvidos, baseados nas mais diferentes premissas, com o intuito de reproduzir numericamente o fenômeno físico como um todo. Tradicionalmente, mé- todos numéricos para reprodução de fenômenos físicos, como a geração de ondas, baseiam-se na discretização do domínio de cálculo em malhas. Entretanto, frente à limitação de representa- ção de domínios de geometrias complexas, descontinuidades ou superfícies móveis, aplicam-se algumas técnicas para contornar tais inconvenientes, como a reestruturação de malha ou ma- lhas adaptativas. Dessa forma, é possível tratar as descontinuidades do domínio, mas que, em contrapartida, aumentam consideravelmente o tempo de processamento além de atrelar a confi- abilidade da solução obtida ao algoritmo de reestruturação de malha. Novos métodos numéricos, desvencilhados do paradigma de discretização da malha, vêm ganhando força nos últimos anos, principalmente pela capacidade em representar fenôme- nos altamente não-lineares e a naturalidade com a qual trata descontinuidades. Estes métodos são conhecidos como métodos de partículas. Os métodos de partículas podem ser definidos, de uma maneira geral, como métodos numéricos de solução de um problema físico que não necessitam de malhas no domínio a ser simulado (meshfree). Este domínio é representado por um conjunto de partículas, cuja dinâmica é determinada de acordo com as equações de balanço (massa e quantidade de movimento), 1.1. OBJETIVO 29 escritas na forma Lagrangiana. Dentre os métodos de partículas destacam-se o SPH (Smoothed Particle Hydrodymanics) e o MPS (Moving Particle Semi-implicit Method). Frente à possibilidade de eliminar as restrições decorrentes da utilização de malhas em domínios de geometrias cada vez mais complexas e ao desenvolvimento recente, eventos e pu- blicações internacionais acerca dos métodos de partículas vêm aumentando significativamente. Esse fato reforça e elucida o interesse da comunidade científica internacional na aplicação dos métodos de partículas aos problemas de Engenharia e sua potencialidade. No entanto, essa ten- dência internacional, no tocante a publicações e formação de recursos humanos, é ainda muito incipiente no cenário nacional. 1.1 OBJETIVO Diante do exposto, pretende-se lançar mão de um código próprio, desenvolvido e testado no âmbito desta Tese de Doutoramento, com base no método SPH, para simular fenômenos que envolvam geração de ondas, com especial atenção aos problemas de ondas de impacto. O código desenvolvido será confrontado com resultados experimentais, analíticos e numéricos, quando houver disposição, de modo a aferir seu desempenho global. 1.2 ORGANIZAÇÃO DO TRABALHO Este documento está organizado em 8 capítulos e duas partes. A primeira parte é dedi- cada à introdução, que compreende a revisão bibliográfica e a descrição dos métodos numéricos. A segunda parte trata especificamente dos estudos de caso efetuados, ou seja, das simulações numéricas de problemas de Engenharia. É dada uma breve descrição do que se espera, em cada capítulo, na sequência. O capítulo 2 é dedicado ao Estado da Arte, em que será apresentado os tipos de mé- todos numéricos, culminando com os métodos particulados. São observados os dois métodos particulados cogitados para serem utilizados no âmbito deste trabalho, sendo o primeiro deles, o MPS, detalhado no capítulo 3. Na sequência, no capítulo 4 é mostrado em detalhes o mo- delo Lagrangiano SPH. Por se tratar de uma técnica relativamente recente, foi incluído todo o equacionamento utilizado e o passo a passo necessário ao seu completo entendimento. Já na Parte II, no capítulo 5, ilustram-se algumas simulações de casos teste, normal- mente análises de desempenho e validação do código numérico. Utilizam-se, para tal fim, a simulação do problema de ruptura de barragens (seção 5.1, p. 59), a geração, propagação e quebra de ondas (seção 5.2, p. 62) e o esvaziamento de um reservatório (seção 5.3, p. 67). O escoamento é considerado, no capítulo 5, bidimensional e ideal, regido pela equação de Euler. Os escoamentos reais são abordados no capítulo 6, através de estudos de caso clássicos em escoamentos de reologia newtoniana, como o Poiseuille em superfície livre (seção 6.2.2, p. 78), o Poiseuille plano (seção 6.2.3, p. 80) e o escoamento Couette (seção 6.2.4, p. 81). 30 Capítulo 1. INTRODUÇÃO Análises transientes e estacionárias são realizadas, onde é possível observar o comportamento do código SPH implementado. Posteriormente, é apresentado o equacionamento do escoamento não newtoniano, a problemática e estudos de caso. No capítulo 7 aborda-se a geração de ondas devido ao impacto de massas sólidas em meio líquido em repouso. É feita uma breve revisão do problema e simulações com material indeformável (seção 7.2.5.1, p. 96) e fragmentado (seção 7.2.5.3, p. 102), comparando com resultados da literatura. As discussões e conclusões a respeito do trabalho em geral são apresentadas no capítulo 8. Faz-se um adendo posterior com as perspectivas para o encaminhamento do trabalho e possí- veis mudanças ou melhorias são propostas, com objetivo de continuação da linha de pesquisa. São referenciados, também, os trabalhos publicados pela equipe de trabalho, com enfoque na temática abordada. Em seguida, apresentam-se a bibliografia citada no texto e apêndices. 31 2 REVISÃO BIBLIOGRÁFICA Dos vários fenômenos físicos observados pelo homem, se é capaz de modelar apenas a parte minoritária, sendo traduzida em equações matemáticas. Estabelecidas as equações mate- máticas e hipóteses, segue-se a tentativa de resolvê-las. Muitas vezes, a solução analítica não é possível para uma determinada situação, em decorrência de dificuldades diversas (topologia complexa, não-linearidades, singularidades, entre outras). Dessa forma, recorre-se a métodos de solução alternativos, como a modelagem numérica. Na modelagem numérica, aplicam-se alguns conceitos às equações originais (diferen- ciais parciais) que regem o fenômeno físico, transformando-as em equações algébricas, para que possam ser resolvidas de forma incremental, em passos. Vistos dessa forma, pode-se dizer que os métodos numéricos nada mais são do que uma ferramenta para a solução de equações algébricas. Existem, no entanto, várias categorias de métodos numéricos, como será visto na sequência. 2.1 FORMULAÇÃO EULERIANA × LAGRANGIANA Existem duas maneiras diferentes de descrever um sistema físico: a Euleriana e a La- grangiana. Essas duas visões distintas são importantes, uma vez que influenciam diretamente no sistema de equações resultante de um dado problema. A formulação Lagrangiana acompanha o movimento de um elemento, espacial e tem- poralmente. O comportamento global do sistema é conhecido quando é determinada a trajetória de todos os elementos do domínio. Já na formulação Euleriana observa-se a variação das gran- dezas pontualmente (na verdade, em uma pequena região do domínio). Os elementos não são acompanhados e o comportamento global do sistema é determinado quando o campo estudado é conhecido em todo o domínio. Existe também uma aproximação mista, chamada ALE, que figura como alternativa para algoritmos de malhas adaptativas. 2.2 MÉTODOS NUMÉRICOS TRADICIONAIS Chamam-se, neste trabalho, de métodos numéricos tradicionais aqueles que necessitam de malha e que, geralmente, seguem a formulação Euleriana. Em termos de métodos numéricos de resolução do problema, atualmente, pode-se apontar como os mais comuns: o Método das Diferenças Finitas (MDF), Volumes Finitos (MVF), Elementos Finitos (MEF), Elementos de Contorno (MEC). A utilização do MDF é adequada para grandes domínios e com requerimentos interme- 32 Capítulo 2. REVISÃO BIBLIOGRÁFICA diários de memória. Entretanto, para compensar a falta de precisão, o domínio deve ser dividido em pequenos elementos, resultando em uma malha densa e estruturada. O MVF pode ser consi- derado como uma variação do MDF, sendo a principal diferença o balanço integral da grandeza avaliada em cada célula do domínio numérico, além de poder ser escrito na forma conservativa e ser mais flexível com o uso de malhas adaptativas. Cabe ainda lembrar o VOF (Volume of fluid), técnica muito utilizada na dinâmica dos fluidos computacional para determinar a posição da superfície livre. Kleefsman et al. (2005), lançam mão da técnica VOF associada ao MVF para determinar a superfície livre, regida pela equação de Navier-Stokes, para simular problemas de impacto em meio liquído, como o embarque de água no convés (green water), analisado como um problema tipo ruptura de barragem com presença de obstáculo posterior, e o slamming de um diedro e cilindro rígidos incidindo sobre meio líquido em repouso. Cheng e Arai (2002) ana- lisam o efeito do sloshing tridimensionalmente em seu modelo numérico em diferenças finitas juntamente com a técnica Surf para determinação da posição da superfície livre. Cheng (1995) modela numericamente em duas dimensões o slamming utilizando as equações de Euler (em coordenadas generalizadas), considerando a separação do escoamento do corpo incidente e ação gravitacional. As equações são discretizadas através de um esquema em diferenças finitas de forma explícita. Maciel (1991) analisa a queda de elementos sólidos in- deformáveis e fragmentados, analisando experimental, numérica e analiticamente o fenômeno de avalanches nos grandes lagos europeus. O autor reproduz a energia adquirida pelo fluido em termos da altura da onda solitária gerada, analisando as teorias de Saint-Venant, Boussinesq e Serre (MACIEL, 1991). Nascimento (2001) usa um modelo numérico (MVF) com base nas equações de Serre para estimar a energia da onda formada pela incidência de um bloco, conside- rado indeformável, que desliza em um plano inclinado e incide num meio líquido em repouso. Mais recentemente, Souza (2007) perfaz o mesmo caminho, agora com material fragmentado (granular) e equipamentos de elevada precisão tanto na determinação do campo orbital quanto nas alturas de onda. Um estudo analítico do impacto hidrodinâmico em águas restritas é rea- lizado por Korobin (1999), através de métodos assintóticos, analisando um objeto tipo caixa. Oliver (2002) também estuda o impacto hidrodinâmico em águas rasas de forma algébrica, res- saltando a importância em conhecer a mecânica desse fenômeno complexo. O MEF possui boa acurácia com poucos elementos, entretanto o método tem a desvan- tagem de alto consumo de memória. O MEF é usado normalmente para descrever o comporta- mento de sólidos. Mitra e Sinhamahapatra (2005) usam o MEF para descrever o comportamento da superfície livre submetida ao sloshing em um tanque retangular usando uma formulação em termos da pressão. Bereznitski (2003) utiliza a formulação MEF tanto para representar a es- trutura quanto o fluido em seu estudo sobre a influência da hidroelasticidade nos impactos hidrodinâmicos. O MEC tem a vantagem de diminuir a dimensão do problema, pois precisa apenas en- contrar a solução no contorno do domínio. Tanizawa (1998) utiliza a formulação MEC no do- 2.3. DIFERENÇAS ENTRE MÉTODOS TRADICIONAIS E MÉTODOS SEM MALHA (MSM) 33 mínio fluido em seu modelo totalmente não-linear, com equações que se baseiam na velocidade e aceleração, em problemas de impacto hidrodinâmico em superfícies elásticas. 2.3 DIFERENÇAS ENTRE MÉTODOS TRADICIONAIS E MÉTODOS SEM MALHA (M- SM) É evidente que, em relação aos métodos numéricos mais tradicionais, como o MDF e MEF, a diferença fundamental dos MSM é que não há necessidade da criação de uma malha. Em geral, nos MSM, o domínio físico é definido por uma porção de partículas, podendo ser geradas de forma aleatória e que não possuem qualquer tipo de conectividade. Entretanto, alguns MSM como o Element-free Garlekin Method (EFGM) se apoiam em células para fazer a integração do sistema matricial (BELYTSCHKO; LU; GU, 1994 apud LIU, 2003), não se caracterizando, portanto, essencialmente como um método de partícula. Há ainda de se ressaltar que, na criação de malhas, as geometrias curvas, por exemplo, são representadas por elementos compostos por segmentos de reta com topologia bem definida e que acompanham a fronteira. Portanto, a equivalência entre o domínio físico e o discretizado é maior quando há mais elementos. Nos MSM, a fronteira pode ser representada por uma série de partículas, sem conectividade entre elas. Pode-se perceber, deste fato, de que a forma da mais complexa geometria pode ser representada sem maiores dificuldades. Na figura 1 é ilustrado um problema dificilmente tratado pelos métodos numéricos tradicionais, com variação brusca da superfície livre, um sólido com movimento livre (três graus de liberdade no plano) e uma fronteira em arco de circunferência. As condições de contorno são fundamentais para que seja obtida a solução do problema proposto. Nos métodos tradicionais, há uma equivalência entre a condição de contorno mate- mática e numérica, sendo uma questão de definir o comportamento das variáveis nas fronteiras do domínio computacional. O mesmo artifício não é tão direto em MSM, até porque não há partículas fora do domínio, o que resulta em diversas técnicas para lidar com essas dificuldades. No caso particular do SPH, os maiores exemplos de condições de contorno são as fronteiras reativas (boundary forces) e as partículas fantasma (ghost particle). Utilizando as condições de contorno do tipo boundary forces, especial atenção deve ser dada às condições iniciais de si- mulação. Em outras palavras, deve-se garantir, nos instantes iniciais, o equilíbrio das partículas. Isso significa que a posição inicial das partículas tem influência no resultado final, diferente do que foi preconizado inicialmente (MONAGHAN, 2006). Condições iniciais também afetam as simulações com partículas fantasma (LACHAMP, 2003). O método SPH em particular aproxima um campo qualquer em um domínio definido por uma mescla de pontos de discretização (partículas). Para cada uma dessas partículas, pode- se escrever uma formulação integral do campo. Entretanto, quando esses métodos são utilizados para resolver problemas de impacto hidrodinâmico, em geral, os esforços calculados apresen- tam oscilações. Especula-se que a origem de tais oscilações possam decorrer dos erros maiores 34 Capítulo 2. REVISÃO BIBLIOGRÁFICA Figura 1 – Exemplo de problema complexo: uma comparação entre um experimento (à es- querda) e o modelo numérico SPH desenvolvido pelo autor (à direita), em instantes aproxima- damente iguais. Fonte: Do próprio autor. 2.4. PASSOS GERAIS PARA RESOLVER UM PROBLEMA USANDO MSM 35 que o SPH comete na aproximação integral de derivadas, em relação à aproximação integral de funções. Dessa forma, o cálculo do gradiente da pressão na equação do balanço de quantidade de movimento incorre em um erro naturalmente maior, que aliado a uma distribuição irregu- lar de massa específica, e por conseguinte pressão, possa resultar nas oscilações de esforços observadas. 2.4 PASSOS GERAIS PARA RESOLVER UM PROBLEMA USANDO MSM Uma vez definidas, de maneira geral, as diferenças e semelhanças entre os métodos tra- dicionais e os MSM, descrevem-se os passos para solução de problemas usando MSM. Primei- ramente, define-se o domínio a ser simulado através da distribuição de partículas. A densidade de partículas pode ser constante ou não, e algumas vezes é necessário que existam mais par- tículas em regiões de interesse. A partícula contém informações diversas como massa, massa específica, velocidade, posição, energia interna, entre outras, dependendo do problema. A relação entre as diversas partículas se dá por meio de uma função peso, que apresenta certas propriedades interessantes e que depende essencialmente do inverso da distância entre as partículas. Desse modo, a maior influência no cálculo das grandezas envolvidas está condicio- nada às partículas mais próximas. Por esse motivo, geralmente limita-se o raio de atuação da função peso a uma distância conhecida da partícula em questão. Essa distância pode ser cha- mada de domínio de suporte (rs, support domain), e deve ser escolhida de modo a garantir uma área adequada para interpolação. Normalmente, o domínio de suporte é escolhido com base na distância característica da simulação, que pode ser tomada como uma fração do espaçamento médio entre partículas (Δx), ou seja, rs = κ Δx. Algumas das propriedades que as funções peso apresentam são comentadas na sequên- cia. A mais importante e comum aos MSM é em relação à unicidade (equação 1), ou seja: ∑ j Wj(�x) = 1. (1) Outras propriedades são desejáveis, e variam de método para método. Uma delas é uma similaridade entre a função peso e a função Delta de Dirac, por exemplo. Uma vez discretizadas, forma-se o sistema de equações a ser resolvido. Normalmente, para problemas de dinâmica dos fluidos, o sistema de equações resultante é composto por equa- ções diferenciais em relação ao tempo, que podem ser resolvidas usando métodos dependendo a abordagem (explícita ou implícita). Não é raro se deparar (e principalmente quando lida-se com problemas implícitos) com o sistema linear do tipo Ax = b, onde a incógnita x é a entrada para o passo de tempo seguinte. Neste caso, há métodos específicos para a solução do sistema linear, que podem ser iterativos ou diretos. Exemplos de métodos diretos são a eliminação de Gauss e a decomposição matricial 36 Capítulo 2. REVISÃO BIBLIOGRÁFICA LU. Já como métodos iterativos, pode-se citar o método de Gauss e suas variantes (Gauss- Jacobi e Gauss-Seidel), Successive over-relaxation method, entre outros. Independentemente da abordagem explícita ou implícita, com a solução do sistema de equações discretizado, obtém-se a variação temporal dos deslocamentos, velocidades e acelerações das partículas. 2.5 PRINCÍPIOS DOS DIFERENTES MÉTODOS NUMÉRICOS O desenvolvimento de um método numérico passa pelo estabelecimento de um princípio matemático e de certas hipóteses de simplificação. Tomando-se inicialmente o MDF: nesse mé- todo, as equações parciais são aproximadas por equações algébricas embasadas no conceito de derivação, através da expansão da função em séries de Taylor. O MEF, por exemplo, baseia-se tradicionalmente nos princípios dos trabalhos virtuais. Da mesma forma, os MSM baseiam-se em diferentes princípios para sua construção. Pode-se dizer também que, quando um método numérico baseia-se nos métodos residuais, por exemplo, suas equações estão escritas na for- mulação fraca (weak form). Os métodos baseados em expansões de Taylor são baseados na formulação forte (strong form). O ideal, em termos de modelagem, seria a obtenção da solução das equações escritas segundo a formulação forte. Esse fato, no entanto, não é possível em alguns problemas de Engenharia. Por esse motivo, uma das alternativas é alterar as equações de balanço do problema em questão, fazendo com que elas não sejam mais totalmente consistentes. Dessa forma, obtém- se soluções fracas em relação a funções ou vetores teste, o que resulta na formulação fraca. Em relação aos MSM, existem aqueles escritos na formulação fraca e forte. Para desen- volver o trabalho de doutoramento em tela, várias alternativas de MSM foram cogitadas para utilização como ferramenta na análise de problemas de geração de ondas. Era desejável que o método adotado seguisse a formulação forte, como por exemplo o MPS e o SPH. Embora o SPH tenha sido criado para resolver problemas no domínio astrofísico (LUCY, 1977; GINGOLD; MONAGHAN, 1977), este método tem sido aplicado em diversas áreas, no- tadamente a mecânica dos fluidos (COLAGROSSI, 2004; LACHAMP, 2003; MONAGHAN, 1982; MORRIS, 1996; OGER et al., 2006; VILA, 1998). Já o MPS foi desenvolvido para li- dar especificamente com escoamentos incompressíveis (KOSHIZUKA; NOBE; OKA, 1998), sendo bastante difundido na Engenharia Naval. 37 3 O MÉTODO MPS O método MPS foi concebido originalmente para lidar com escoamentos incompressí- veis. Em termos gerais, pode-se dizer que o problema da geração de ondas por impacto hidrodi- nâmico apresenta baixa dependência com a compressibilidade do fluido, o que tira ainda maior proveito da arquitetura desse método. Sendo assim, as equações motrizes utilizadas no MPS são as equações da continuidade e de Navier-Stokes, com a hipótese de escoamento incompressível: �∇ ·�u = 0, (2) d�u dt =− 1 ρ �∇p+ν�∇2�u+�g. (3) Adota-se, para restringir as interações entre as partículas a uma região de raio finito (rs), uma função peso ϖ , definida como: ϖ(r) = { rs ri j −1 ri j < rs, 0 ri j ≥ rs. (4) A função peso pode, no entanto, apresentar formulações diferentes. Um estudo nesse sentido foi empreendido por Ataie-Ashtiani e Farhadi (2006), comparando os resultados experimentais e numéricos resultantes de seis formulações diferentes para a função peso, avaliado para o problema tipo ruptura de barragens. Um indicativo da massa específica da partícula, Ni, é definido por: Ni = ∑ j =i ϖ(r ji). (5) Pela definição do gradiente, pode-se perceber que o operador nabla aplicado a um campo escalar φ qualquer entre duas partículas i e j, é: �∇φ = φ j −φi r2 ji (�r j −�ri). (6) O vetor gradiente da partícula i é obtido medianizando os vetores gradiente segundo a função peso: �∇φi = D Ni ∑ j =i [ (φ j −φi)(�r j −�ri) r2 ji ϖ(r ji) ] . (7) Já o operador laplaciano é dado por: ∇2φi = 2D λ N0 ∑ j =i (φ j −φi)ϖ(r ji), (8) 38 Capítulo 3. O MÉTODO MPS com: λ = ∫ v ϖ(r)r2 dv∫ v ϖ(r)dv . (9) Para resolver o sistema de equações governantes, utiliza-se um algoritmo semi-implícito. Primeiramente, no passo explícito, a dinâmica das partículas é determinada pela solução da equação de conservação de quantidade de movimento: (�u)∗i = (�u)t i +Δt (ν∇2(�u)t i +�gi), (10) (�r)∗i = (�r)t i +Δt (�u)∗i . (11) No passo seguinte, a equação da continuidade é resolvida implicitamente: �∇ · (�u),i = − 1 Δt ρ t+Δt i −ρ∗ i ρ0 =− 1 Δt Nt+Δt i −N∗ i N0 , (12) (�u),i = −Δt ρ (�∇pt+Δt i ). (13) com (�u),i = (�u)t+Δt i +(�u)∗i . Na equação 12, substitui-se a massa específica pelo seu respectivo indicativo (equação 5). Combinando as duas equações acima resulta: (∇2 p)t+Δt i = ρ Δt2 [ Nt+Δt i −N∗ i N0 ] . (14) Adotando que o indicativo de massa específica do fluido é uma função linear da massa específica, tem-se: Nt+Δt i N0 = ρ t+Δt i ρ0 , (15) desconsiderando os efeitos de compressibilidade do fluido. Substituindo a equação 15 na equa- ção 14 chega-se a equação de pressão final que é resolvida pelo método: (∇2 p)t+Δt i =− ρ Δt2 [ N∗ i −N0 N0 ] . (16) 39 4 O MÉTODO SPH O método SPH foi criado inicialmente para resolver problemas astrofísicos no espaço tridimensional (LUCY, 1977; GINGOLD; MONAGHAN, 1977). No entanto, a simplicidade com que fenômenos complexos são modelados faz do SPH uma ferramenta de solução inte- ressante, tendo sido estudada extensivamente e estendida para diversos problemas e áreas de atuação (VILA, 1998). O SPH, como é um método sem malha, lagrangiano e de partícula, possui características particulares. Sem dúvida, possui vantagens especiais sobre os métodos numéricos tradicionais, baseados em malhas, sendo a mais significativa a natureza adaptativa do SPH. Essa adaptabili- dade é alcançada em um estágio inicial na aproximação da variável do campo, que é realizada em cada passo de tempo baseada na quantidade atual e local de partículas. O significado da sigla SPH exprime a filosofia do método. O primeiro termo Smoothed representa a suavização ou medianização segundo o peso das variáveis do campo na vizinhança da partícula, visando principalmente a estabilidade do método. O termo Particle refere-se à utilização de partículas para discretização do domínio físico. O terceiro termo Hydrodynamics representa o nicho de atuação do método. Assim, a combinação da adaptabilidade, natureza particulada e lagrangiana leva o método SPH a ser aplicado em diferentes áreas da engenharia e ciência. De certa maneira, o termo Hydrodynamics pode ser interpretado como mecânica de forma geral. Em algumas literaturas (KUM; HOOVER; POSCH, 1995 apud LIU, 2003) esse método é chamado de Smoothed Particle Mechanics. 4.1 FUNDAMENTOS DO MÉTODO SPH Do ponto de vista computacional, representa-se o fluido por uma porção de partículas evoluindo com a velocidade do escoamento. Cada partícula representa um ponto de interpolação na qual todas as propriedades do fluido são conhecidas. Exprimindo então a função φ como um campo variável qualquer (escalar, vetorial ou tensorial), a seguinte igualdade verifica: φ(�r) = ∫ φ(�r j)δ (�r−�r j)d�r j (17) Se a função delta de Dirac for substituída por uma função W (�r−�r j,h), que satisfaça as seguintes condições: ∫ W (‖�r−�r j‖,h)d�r j = 1 (18) lim h→0 W (‖�r−�r j‖,h) = δ (�r−�r j) (19) 40 Capítulo 4. O MÉTODO SPH A função φ pode, então, ser reescrita como: φ(�r)≈ ∫ φ(�r j)W (‖�r−�r j‖,h)d�r j (20) Na forma apresentada pela equação 20, a função φ(�r) pode ser considerada como sendo uma regularização da função φ original. O termo W inserido, cujas propriedades assemelham- se às do delta de Dirac, justificam a nomenclatura núcleo de suavização (smoothing kernel) e o parâmetro h é o domínio de suporte, que no SPH recebe a denominação de comprimento de suavização (smoothing lenght) do núcleo. Aplicando a equação 20 no domínio discretizado, tem-se: 〈φ(�r)〉= ∑ j φ j m j ρ j W (‖�r−�r j‖,h) (21) no qual, de maneira geral, o valor de φ em�r j (ou na partícula j) é denotado por φ j. Pode-se ainda aplicar o gradiente à função φ , utilizando novamente a equação 20: �∇φ(�r)≈ ∫ �∇φ(�r j)W (‖�r−�r j‖,h)d�r j (22) Integrando a equação 22 por partes, desprezando os termos de superfície e aplicando as propriedades do núcleo W (OGER et al., 2007) que serão discutidas posteriormente (na seção 4.2, p. 41), resulta: �∇φ(�r)≈ ∫ φ(�r j)�∇W (‖�r−�r j‖,h)d�r j (23) Que pode ser escrito como: 〈�∇φ(�r)〉= ∑ j φ j m j ρ j �∇W (‖�r−�r j‖,h) (24) No entanto, utiliza-se comumente a equação 25 (MONAGHAN, 2005): ρ�∇φ(�r) = �∇(ρ φ(�r))−φ(�r)�∇ρ (25) Que dá origem a uma das expressões para o gradiente de φ que pode ser utilizada no SPH: 〈�∇φ(�r)〉i = 1 ρi ∑ j m j ( φ j −φi ) �∇W (ri j,h) (26) Da mesma forma, para o divergente do campo vetorial�u, pode-se obter diretamente: 〈�∇ ·�u〉i = ∑ j m j�u j ·�∇W (ri j,h) (27) Novamente, opta-se por um rearranjo no divergente, escrevendo-o como (MONAGHAN, 2005): �∇ ·�u = 1 ρ [�∇ · (ρ�u)−�u ·�∇ρ] (28) 4.2. A FUNÇÃO NÚCLEO DE SUAVIZAÇÃO W 41 Que resulta na proposta de aproximação do divergente de uma função vetorial �u, muito comum na formulação SPH: 〈�∇ ·�u〉i = 1 ρi ∑ j m j (�u j −�ui) ·�∇iW (ri j,h) (29) sendo que o termo �∇iW (ri j,h) representa o gradiente de W (ri j,h) em relação às coordenadas da partícula i. As representações do divergente e gradiente de uma função qualquer serão utilizadas para reescrever as equações de balanço (continuidade e quantidade de movimento) segundo a forma SPH. De modo a facilitar, emprega-se a notação: Wi j = W (ri j,h), φi j = φi − φ j e φ̄i j = (φi+φ j)/2. Os símbolos indicando funções aproximadas (〈·〉) também serão omitidos, de agora em diante, considerando-os implícitos nas equações. 4.2 A FUNÇÃO NÚCLEO DE SUAVIZAÇÃO W O núcleo de suavização é uma função analítica e contínua e tem a responsabilidade de interpolar uma grandeza qualquer entre partículas, limitada pela máxima distância K h, K con- tante. Pode-se dizer que o núcleo de suavização tem o mesmo papel dos esquemas de discreti- zação nos métodos tipo MDF (central, upwind, etc.). Sendo respeitadas as restrições impostas pelas equações 18 e 19, qualquer função pode ser candidata a representar W . A regulariza- ção gaussiana (equação 30) pode ser considerada a primeira categoria de núcleo de suavização aplicada (MONAGHAN; KOS; ISSA, 2003): W G ( r h ) = 1√ π exp [ − ( r h )2 ] (30) Para facilitar a notação, comumente adota-se s = r/h. Existem vantagens na adoção do núcleo gaussiano, sendo uma delas obtida pelo com- portamento da função exponencial, no sentido que derivadas do núcleo pode ser escrita em função do próprio núcleo: dW G ds =−2sW G (31) Além disso, Morris (1996) aponta que a transformada de Fourier do núcleo gaussiano é gaussiano e que há boa estabilidade do método como um todo quando o núcleo gaussiano é utilizado. No entanto, a falta de suporte compacto faz com que todas as partículas do domínio influenciem umas as outras (ou seja, K é infinito), mesmo que as contribuições sejam pequenas para partículas afastadas. Um outro tipo de núcleo bastante utilizado são as funções splines, sendo a cúbica bidi- 42 Capítulo 4. O MÉTODO SPH mensional dada pela equação 32: W SC(s) = 10 7πh2 ⎧⎪⎪⎨ ⎪⎪⎩ 2 3 − s2 + s3 2 , se 0 < s ≤ 1 (2−s)3 6 , se 1 < s ≤ 2 0, se s > 2 (32) Uma das maiores vantagens na utilização do núcleo regido pela spline cúbica é que a função é exatamente nula para s > 2, ou seja, para uma distância superior a 2h, como pode ser visto na equação 32. Essa propriedade é o suporte compacto, que reduz drasticamente o tempo computacional, já que a contribuição das partículas é limitada a um raio 2h. No entanto, em algumas aplicações, pode-se notar acentuados efeitos dispersivos (MORRIS, 1996). A figura 2 ilustra a aplicação do núcleo de suavização, no caso em que o raio de interação é limitado em 2h. Figura 2 – À esquerda, o raio de ação da função W e à direita os valores da grandeza inter- polados, representados pelas barras verticais cuja altura é proporcional ao inverso da distância ri j. Fonte: Do próprio autor. Independente da formulação do núcleo de suavização, é importante lembrar que a trans- formação do domínio contínuo para o discreto passa por estabelecer as distâncias entre as par- tículas intervenientes, uma vez que W →W (ri j,h). O estabelecimento das partículas vizinhas, ou simplesmente da vizinhança, exige um esforço computacional, que pode ser reduzido pela propriedade de suporte compacto de alguns núcleos de suavização. Entretanto, para que seja ex- plorada as vantagens do suporte compacto, algumas alterações são exigidas no SPH. Em outras palavras, não há vantagens significativas em usar o núcleo regido por splines se é necessário varrer todas as partículas para verificar se a distância até a partícula considerada é menor que 2h. Para contornar esse inconveniente, Morris (1996) apresenta um algoritmo de busca de vizi- nhos, baseado na subdivisão do domínio em células quadradas, de lado 2h. Conhecendo-se, de 4.3. COMPRIMENTO DE SUAVIZAÇÃO 43 antemão, quais partículas pertencem a cada célula (através de uma lista de ligação), é possível restringir a pesquisa das partículas vizinhas apenas às células adjacentes (células pintadas na figura 3). Esse artifício reduz drasticamente o tempo computacional exigido para a construção das vizinhanças das partículas. Há outras técnicas, como a tabela Verlet (COLAGROSSI, 2004, Verlet table), que consiste na introdução de um segundo raio, maior que o original apenas por um valor Δ. Garante-se que, para um dado número de passos de tempo, uma partícula, mesmo com a máxima velocidade permitida, não ultrapassaria a região ao qual pertence. Dessa forma, otimiza-se ainda mais a busca por vizinhos, garantindo que durante os números de passo de tempo estabelecidos, todos os vizinhos da partícula em questão estarão incluídos. Figura 3 – Subdivisão do domínio computacional para diminuir o número de partículas de busca na criação da lista de vizinhança. À esquerda, vê-se todo o domínio computacional, en- quanto à direita vê-se o detalhe da região de interesse. Fonte: Do próprio autor. A diferenciabilidade do núcleo de suavização é um fator importante, uma vez que a aproximação do gradiente de uma função qualquer depende exclusivamente das derivadas da função de suavização, como mostram as equações 26 e 29. Há um consenso geral de que o erro no cálculo de derivadas de funções (equação 26) é maior do que a aproximação de funções (equação 21) (MONAGHAN, 2005). A figura 4 ilustra os núcleos gaussiano e spline cúbico. 4.3 COMPRIMENTO DE SUAVIZAÇÃO O comprimento de suavização h está intimamente relacionado à resolução numérica desejada, e é uma denominação particular do SPH para o domínio de suporte (rs). Se o escoa- mento simulado não apresenta regiões com grande variação de massa específica, o comprimento h pode-se manter constante. Contudo, para simular choques ou impactos, onde há variação sú- bita da massa específica localmente, é preciso variar h de acordo, de tal forma que o número de partículas vizinhas se mantenha próximo de um valor constante. 44 Capítulo 4. O MÉTODO SPH Figura 4 – Função de suavização ou núcleo W . Fonte: Do próprio autor. Toma-se h como uma proporção do espaçamento da partículas inicial (Δx). Normal- mente, em problemas bidimensionais, h = 1,2Δx. Para esse valor de h, e considerando o su- porte compacto do núcleo de raio 2h, espera-se a presença de, em média, 20 partículas vizinhas (MORRIS, 1996). Propõe-se, como aproximação inicial, uma expressão que relacione diretamente h com a variação da massa específica (equação 33). Essa expressão, todavia, necessita estabelecer previamente a massa específica da partícula i, cujo cálculo depende de hi. hi = h0 ( ρ0 ρi ) 1 d (33) Desse modo, para estabelecer uma expressão para h, deriva-se a equação 33 com relação a t, injetando a equação da continuidade no termo resultante, obtém-se: dh dt = 1 d h�∇ ·�u (34) 4.4 EQUAÇÃO DA QUANTIDADE DE MOVIMENTO A equação da quantidade de movimento: d�ui dt =− ( �∇p ) i ρi (35) 4.5. DESLOCAMENTO DE PARTÍCULAS 45 implica na representação do gradiente de um escalar na formulação SPH. Aplicando a equação 26, resulta: ( �∇p ) i =− 1 ρi ∑ j m j pi j�∇iWi j (36) A expressão 36, contudo, não conserva a quantidade de movimento linear e angular exatamente. Para solucionar esse problema, reescreve-se o gradiente de outra forma, de acordo com a equação 37: �∇p ρ = �∇ ( p ρ ) + p ρ2 �∇ρ (37) A equação resultante para a quantidade de movimento fica: d�ui dt =−∑ j m j ( pi + p j ρi ρ j ) �∇iWi j (38) que conserva exatamente a quantidade de movimento linear e angular. No entanto, existem ou- tras maneiras de reescrever o gradiente que conservam a quantidade de movimento exatamente (LACHAMP, 2003; MORRIS, 1996). 4.5 DESLOCAMENTO DE PARTÍCULAS As partículas são deslocadas de acordo com: d�ri dt =�ui (39) Existem métodos que melhoram a estabilidade da equação 39. Um desses métodos é o XSPH, que adiciona o seguinte termo à equação 39: Δ�ui = ε ∑ j m j �ui j ρ̄i j Wi j (40) sendo ε uma constante que varia entre 0 e 1 (valor usual: ε = 0,5). Segundo Monaghan e Kos (1999), essa correção mantém as partículas mais ordenadas, além de prevenir penetração em casos de escoamentos de alta velocidade. Lachamp (2003) ressalta que o termo corretor sugerido (equação 40) não introduz dissipação viscosa suplementar, aumentando ligeiramente, todavia, a dispersão do sistema. 4.6 EQUAÇÃO DA CONTINUIDADE Devido às características do sistema particulado, pode-se dizer que o balanço de massa é dispensável no SPH, uma vez que a massa das partículas é constante e, numa determinada simulação, não há acréscimo nem redução no número de partículas. No entanto, o escoamento de fluidos como a água, no SPH, é considerado como fracamente compressível. Dessa forma, 46 Capítulo 4. O MÉTODO SPH não se garante que a massa específica seja constante, tampouco �∇ ·�u = 0. Ou seja, resolve-se a equação: dρ dt =−ρ�∇ ·�u (41) Como já visto (seção 4.1, p. 39), pode-se obter a aproximação desejada aplicando dire- tamente o conceito do SPH ao divergente (equação 24), que resulta: ρi = ∑ j m j Wi j. (42) A equação 42 conserva exatamente o volume total, uma vez que o número de partícu- las e massa de cada partícula são constantes ao longo da simulação. Mas, principalmente em problemas com superfície livre, nota-se um decaimento acentuado na massa específica nas pro- ximidades da superfície livre, o que resulta em comportamentos indesejáveis nessa região. Para contornar esses inconvenientes, opta-se por um rearranjo na formulação do divergente (equação 28), que aplicado à equação da continuidade resulta: dρi dt = ρi ∑ j m j ρ j �ui j ·�∇iWi j (43) Utilizando a equação 43, o problema de queda repentina na massa específica é con- tornado pois a massa específica pode ser atribuída a cada partícula, fazendo assim com que a variação dρ/dt dependa da aproximação ou afastamento relativo entre as partículas (termo �ui j não nulo). Podem ser aplicadas correções e/ou suavizações à massa específica de tempos em tempos na simulação, para contornar o fato da equação 43 não conservar o volume total. 4.7 CORREÇÕES APLICADAS AO SPH As imposições ao núcleo de suavização (equações 18 e 19), que devem ser obedecidas para que o núcleo de suavização assuma as propriedades do delta de Dirac num limite quando h → 0, podem ser eventualmente violadas. Isto pode acontecer dada uma certa desorganização no posicionamento das partículas, que somam um número finito distribuídas em um domínio discreto. Em outras palavras, não é garantido que o número de vizinhos de cada partícula perma- neça constante, nem tampouco guardem a mesma distância entre partículas. Colagrossi (2004) mostra as equações que normalmente deveriam, mas não são sempre satisfeitas pelo núcleo de suavização: ∑ j m j Wi j/ρ j = 1 ∑ j m j�∇Wi j/ρ j =�0 ∑ j �r j m j Wi j/ρ j =�ri ∑ j m j�r j ·�∇Wi j/ρ j = I. (44) A seguir, são discutidas duas técnicas diferentes para corrigir alguns erros no SPH: A renormalização e a reinicialização da massa específica. 4.7. CORREÇÕES APLICADAS AO SPH 47 4.7.1 Reinicialização da massa específica A reinicialização da massa específica é uma correção aplicada em intervalos de tempo pré-definidos (LAIGLE; LACHAMP; NAAIM, 2007) cujo objetivo é amenizar os efeitos nega- tivos na utilização da equação 43 para a equação da continuidade. A correção envolve aplicar, para todas as partículas, a equação 42, que conserva totalmente o volume na simulação. Entretanto, caso o núcleo de suavização Wi j não seja normalizado, ou seja, ∑ j Wi j = 1, a aplicação da correção insere ainda mais erros, causando a degradação do campo de pressões. Sendo assim, altera-se o núcleo de suavização original, através dos interpoladores de Shepard ou o MLS (Moving Least Squares). 4.7.1.1 Interpolador de Shepard O interpolador de Shepard (SHEPARD, 1968) é uma técnica simples, de baixo consumo de tempo computacional, que permite uma correta interpolação de funções constantes, ou seja, pode-se dizer que é uma técnica de ordem zero. A correção aplicada ao núcleo de suavização é dada por: W IS i j = Wi j ∑l ml Wil/ρl . (45) 4.7.1.2 Interpolador MLS O interpolador MLS apresenta ordem n, sendo capaz de interpolar exatamente um po- linômio de mesma ordem. No entanto, por causa do número crescente de operações e tempo computacional, limita-se a técnica à ordem 1. Para o caso supracitado, o núcleo de suavização é dado por: W MLS i j = ( ζ (0)i +ζ (1)i(x j − xi)+ζ (2)i(y j − yi) ) Wi j, (46) sendo que o vetor ζi é calculado, bidimensionalmente, resolvendo um sistema de equações lineares, Ai jζi = b, sendo: Ai j ⎡ ⎢⎣ ζ (0)i ζ (1)i ζ (2)i ⎤ ⎥⎦= ⎡ ⎢⎣ 1 0 0 ⎤ ⎥⎦ , (47) com: Ai j = ∑ j ⎡ ⎢⎣ 1 x j − xi y j − yi x j − xi (x j − xi) 2 (x j − xi)(y j − yi) y j − yi (x j − xi)(y j − yi) (y j − yi) 2 ⎤ ⎥⎦Wi j m j ρ j . (48) É interessante notar que, em casos onde a partícula tem poucos vizinhos, a matriz Ai j pode ser singular. Como esse fato impossibilita a aplicação do interpolador de forma plena, faz-se uma verificação da possibilidade de redução da matriz, para uma submatriz que possua determinante não-nulo. Esta característica de envelopamento de Ai j constitui uma das caracte- rísticas do interpolador MLS, que baseia-se em uma ordem variável (DILTS, 1999). 48 Capítulo 4. O MÉTODO SPH Como a ordem do interpolador escolhida é 1, logicamente que a submatriz resultante dá origem a um método de ordem zero, que é o interpolador de Shepard (ζ (1)i e ζ (2)i iguais a zero na equação 46). 4.7.2 Renormalização Por meio da renormalização, procura-se satisfazer as restrições ao núcleo de suaviza- ção (equação 18), independente da escolha do comprimento de suavização ou da organização das partículas. Bonet e Lok (1999) destacam as vantagens das correções, tanto no núcleo de suavização quanto no gradiente. Sendo assim, uma função qualquer pode ser calculada como: φi = ∑ j φ j m j W IS i j /ρ j. (49) sendo W IS i j dado pela equação 45. Essa correção nada mais é do que a aplicação do interpolador de Shepard. Quanto ao gradiente, tem-se a equação 50: �∇φi = ∑ j φ j m j ∇̃W IS i j /ρ j. (50) sendo que o termo ∇̃W IS i j é obtido corrigindo o gradiente do núcleo de suavização W IS i j , fazendo: ∇̃W IS i j = Bi j∇W IS i j . O cálculo de ∇W IS i j é obtido aplicando o gradiente à equação 45, que resulta em: ∇W IS i j = ( ∇Wi j ∑ j m j Wi j/ρ j )− ( Wi j ∑ j m j ∇Wi j/ρ j ) ( ∑ j m j Wi j/ρ j )2 (51) Finalmente, para determinar ∇̃W IS i j resta estabelecer Bi j. Essa matriz é obtida impondo a condição ∫ (x j − xi)Wi jdr j = 0 ao gradiente do núcleo de suavização, que resulta em: Bi j = ⎡ ⎣∑ j m j ρ j ∂W IS ∂ r x j(xi−x j) ri j ∑ j m j ρ j ∂W IS ∂ r y j(xi−x j) ri j ∑ j m j ρ j ∂W IS ∂ r x j(yi−y j) ri j ∑ j m j ρ j ∂W IS ∂ r y j(yi−y j) ri j ⎤ ⎦ −1 (52) No trabalho de doutoramento em tela, opta-se apenas pela renormalização do gradiente do kernel. Dentre os motivos principais, destaca-se que a correção mista sugerida por Bonet e Lok (1999) (equação 50) pode ser substituída pela aproximação: �∇φi = ∑ j (φ j −φi)m j ∇Wi j/ρ j, (53) que ainda campos constantes serão calculados exatamente. Adotando a equação 53, a matriz de correção Bi j fica: Bi j = ⎡ ⎢⎣ −∑ j m j ρ j ∂W IS ∂ r (xi−x j) 2 ri j ∑ j m j ρ j ∂W IS ∂ r (xi−x j)(y j−yi) ri j ∑ j m j ρ j ∂W IS ∂ r (yi−y j)(x j−xi) ri j −∑ j m j ρ j ∂W IS ∂ r (yi−y j) 2 ri j ⎤ ⎥⎦ −1 (54) 4.8. LEIS DE ESTADO PARA PRESSÃO 49 4.8 LEIS DE ESTADO PARA PRESSÃO A modelagem da pressão é um ponto delicado da técnica SPH aplicada a escoamentos incompressíveis, devido à falta de controle explícito da massa específica local. Como o SPH converge bem para escoamentos compressíveis, aproxima-se o caso incompressível por um escoamento quase-compressível através de uma equação de estado para a pressão. Essa equação de estado é da forma p = p(ρ) e, além de fechar o sistema de equações do SPH, mantém sua característica essencialmente explícita. Outras técnicas, como o MPS (e também o ISPH, ou SPH Incompressível), dispensam a utilização de uma equação de estado pois apostam na solução implícita da equação da pressão (equação de Poisson, ∇2φ +K = 0, K é uma constante), gerando naturalmente um cenário de incompressibilidade. Cabe lembrar que a aproximação para o caso incompressível está intimamente ligada à escolha da celeridade do som da simulação. Isso porque, em escoamentos com baixo número de Mach (Ma), as variações de massa específica são proporcionais a Ma2, sendo Ma dado por: Ma = u/c (55) Então, restringindo-se as velocidades no fluido para que Ma seja da ordem de O(10−1), garante-se que as variações de massa específica sejam na ordem de 1%. Normalmente, a ce- leridade do som utilizada na simulação segue esse procedimento, isso porque a utilização do valor físico da celeridade do som implicaria em um tempo de processamento muito maior, sem no entanto qualquer ganho significativo. Em outras palavras, o valor de c na simulação deve ser suficientemente alto para reduzir as oscilações de ordem numérica (ex: simulação de esco- amentos ideais) e suficientemente baixo para permitir que o modelo tenha um passo de tempo razoável. Assim, via de regra, considera-se um valor c bem menor do que o seu valor real. É comum adotar c = 10V , sendo V a máxima velocidade encontrada na simulação. Existem diversas propostas para a equação de estado. Uma das formulações mais co- muns é aquela dada pela equação 56, sendo a relação entre pressão e massa específica dada por: pi = c2ρ0 7 [( ρi ρ0 )7 −1 ] (56) Laigle, Lachamp e Naaim (2007) utilizam uma expressão bem simples para a pressão (equação 57), e segundo os autores resultados satisfatórios foram obtidos da análise de escoa- mentos não-Newtonianos sobre obstáculos: pi = p0 + c2(ρi −ρ0) (57) 50 Capítulo 4. O MÉTODO SPH 4.9 VISCOSIDADE ARTIFICIAL Observa-se, em problemas de escoamento regidos pela equação de Euler sem qualquer tipo de viscosidade, oscilações que não correspondem à física, geralmente ocasionadas por uma difusão numérica de pequena ordem. Assim como em outros métodos, adiciona-se à equação de quantidade de movimento um termo para anular essas oscilações, correspondendo à conhecida viscosidade numérica. No SPH, vários termos de correção já foram propostos (LACHAMP, 2003; MORRIS, 1996), mas aquele que é mais largamente utilizado é uma espécie de pressão viscosa, dada por: Πi j = ⎧⎨ ⎩ −α c̄i j μ̃i j+β μ̃2 i j ρ̄ , se �ui j ·�ri j < 0 0, se �ui j ·�ri j ≥ 0 (58) com: μ̃i j = h�ui j ·�ri j ri j2 +0,01h2 (59) Introduzindo o termo Πi j (equação 58) na equação de quantidade de movimento (equa- ção 38), resulta: d�ui dt =−∑ j m j ( pi + p j ρi ρ j +Πi j ) �∇iWi j (60) Pode-se perceber pela equação 60 o motivo da denominação de pressão viscosa ou pres- são artificial para o termo Πi j, uma vez que é inserido diretamente no termo de gradiente de pressões da equação de quantidade de movimento. Recentemente, Kajtar e Monaghan (2008) propõem um novo tipo de equação para a viscosidade artificial, que é dada por: Πi j =−α vsig�ui j ·�j ρ̄i jri j (61) Prova-se que as duas formulações para viscosidade artificial (equações 58 e 61) são idênticas (pode ser verificado no capítulo 6, seção 6.1.2, p. 72). Percebe-se, pela maneira como a viscosidade artificial é equacionada, que ela é diferente de zero apenas quando �ui j ·�ri j ≥ 0. Isso ocorre porque quando o produto escalar em destaque é negativo, as partículas estão se aproximando. Nesta situação, a viscosidade artificial é posi- tiva, gerando uma força de repulsão entre as partículas. No caso em que o produto escalar é positivo, as partículas estão se afastando e a viscosidade artificial é negativa, gerando uma força de atração entre as partículas. Dessa forma, é comum considerar a viscosidade artificial apenas quando as partículas aproximam-se umas das outras, anulando seu efeito de atração no caso de distanciamento. 4.10 CONDIÇÕES DE CONTORNO Originalmente, o SPH foi concebido para tratar de problemas no domínio astrofísico, portanto em problemas com condições de fronteira menos restritivas que fronteiras sólidas, por 4.10. CONDIÇÕES DE CONTORNO 51 exemplo. Além disso, não há a necessidade de estabelecer certas condições de contorno. Por exemplo, não é necessário estabelecer a condição dinâmica ou cinemática na superfície livre. Essas condições estão, de certa forma, incorporadas à formulação do método (OGER et al., 2007). Entretanto, na grande maioria dos problemas em Engenharia, lida-se com condições de contorno mais restritivas (impermeabilidade de paredes, condição de aderência, etc). Tais condições de contorno não são satisfeitas automaticamente pelas equações do método. Dessa forma, o tratamento das fronteiras recebe atenção especial e seguem basicamente duas vertentes: a adoção das partículas fantasma (ghost particles) ou fronteiras reativas (boundary forces). 4.10.1 Partículas fantasma No caso das partículas fantasma, quando uma partícula aproxima-se da fronteira, há o espelhamento da partícula, criando uma partícula fictícia com a mesma distância normal à pa- rede. Essa partícula fictícia, ou partícula fantasma, possui as mesmas características da partícula real, entretanto o sentido da velocidade segue a condição de contorno da parede. Esse artifício aumenta a vizinhança da partícula próxima à parede, garantindo unicidade do somatório do nú- cleo de suavização (equação 18). A figura 5 exemplifica como são criadas as partículas fantasma para simular uma condição de contorno. Figura 5 – Partículas fantasma, com condição �ugh ·�t =−�ui ·�t, ∀i. Fonte: Do próprio autor. De acordo com a condição de impenetrabilidade, a componente normal da velocidade da partícula fantasma deve ser contrária à da partícula real. Sendo assim: �ugh ·�n =−�ui ·�n, (62) para fronteiras que não se movem. 52 Capítulo 4. O MÉTODO SPH Já para a componente tangencial da velocidade da partícula fantasma, existem basica- mente três maneiras diferentes de representação, considerando a fronteira em repouso: a) �ugh ·�t =−�ui ·�t, ∀i b) �ugh ·�t = �ui ·�t, ∀i c) �ugh ·�t = 0, ∀i Na primeira forma, a componente tangencial da velocidade da partícula fantasma é igual e oposta à da partícula fluida. Esta seria a melhor tradução da condição de contorno do tipo aderência (no-slip). Já a segunda forma lança mão da simetria, ou espelhamento das proprieda- des. Esta condição de contorno representa bem a condição de deslizamento (free-slip) (COLA- GROSSI, 2004). Finalmente, a terceira forma de representar a componente tangencial da velocidade das partículas fantasma é considerá-la com a mesma velocidade da fronteira (velocidade relativa nula), qualquer que seja i. Embora não existe equivalência matemática para esse conceito, pode- se considerar que o efeito que essa imposição causa ao escoamento é da aderência. Rodriguez-Paz e Bonet (2004) abordam a condição de aderência para fluidos Bingha- mianos, utilizando o SPH, como uma compatibilidade de forças na fronteira, adicionando o atrito do fundo de acordo com o atrito interno do material escoante. Além disso, os autores também discorrem a respeito da possibilidade das partículas escaparem do domínio, passando através de uma fronteira impermeável. Nesta situação, considera-se que há uma condição de reflexão, de maneira que a partícula é reintroduzida no domínio em uma posição que depende da elasticidade da fronteira. Takeda, Miyama e Sekiya (1994) utilizam uma aproximação que mescla a atribuição de velocidades das partículas fantasma com partículas fixas. Os autores, para estabelecer as condições de contorno cinemáticas nas fronteiras do escoamento Poiseuille plano criam, de forma aleatória, partículas fora do domínio fluido, as chamadas partículas imaginárias. Quando uma partícula real tem uma imaginária como vizinha, a velocidade da imaginária é calculada com: �ugh = −dB/dA�ui, onde dB é a distância normal da partícula imaginária gh à fronteira e dA é a distância normal da partícula real i à fronteira. A proposta de Takeda, Miyama e Sekiya (1994) garante uma interpolação linear da velocidade entre as partículas i e gh de tal forma que, na fronteira, a velocidade é nula. Morris, Fox e Zhu (1997) utilizam uma equação semelhante, limitando o valor de �ugh em 1,5�ui. Recentemente, em um estudo feito por Maciá et al. (2011), os autores verificaram o cál- culo do laplaciano da velocidade na fronteira em escoamentos newtonianos no SPH. A despeito dos cálculos incorretos, os autores apontaram que a condição de aderência apresentada por Ta- keda, Miyama e Sekiya (1994) é a que mais se aproxima do valor esperado, dentre as formas usuais de representar esta condição de contorno. 4.10. CONDIÇÕES DE CONTORNO 53 A implementação das partículas fantasma seguiu os preceitos apontados anteriormente, através de uma matriz de transformação. No capítulo 6 serão feitos vários testes numéricos utilizando as partículas fantasma. 4.10.2 Fronteiras reativas No caso de fronteiras reativas, utiliza-se o conceito de atração e repulsão molecular (MONAGHAN; KOS, 1999). Sendo assim, a fronteira é representada por uma linha de partícu- las reais que exercem uma força a qualquer partícula que entrar em seu raio de vizinhança. Em termos de equacionamento, as fronteiras reativas são tratadas como proposto por Peskin (ou, mais precisamente, por Sirovich, já que o autor considera as fronteiras não elásti- cas, diferentemente de Peskin) (MONAGHAN; KAJTAR, 2009). Segundo estes autores, não é feita uma analogia matemática das condição de contorno, mas sim adiciona-se um termo à equação fonte (quantidade de movimento) cujo efeito global seja o mesmo. Em outras palavras, não é feita uma compatibilidade cinemática correspondente à condição de contorno na fronteira (�u normal à parede nula), mas sim uma compatibilidade de forças. Sabe-se que com esse arti- fício, Peskin (2002) vem obtendo sucesso com o Método de Fronteira Imersa na representação de sistemas cardiovasculares e aplicações na biomecânica em geral. Inserindo o termo repre- sentativo da força exercida pela parede na equação de quantidade de movimento (60), tem-se: d�ui dt =−∑ j m j ( pi + p j ρi ρ j +Πi j ) �∇iWi j +∑ k ( �fik −mkΠik �∇iWik ) (63) sendo que a soma em j abrange as partículas fluidas enquanto k as partículas de parede. Ressaltam-se os problemas oriundos do acréscimo de um termo suplementar à equa- ção de quantidade de movimento, em problemas que exigem a estacionaridade como condição inicial. Nessas situações, frequentemente faz-se um aquecimento do modelo, através de um amortecimento artificial, para que as partículas se acomodem e que a posição de equilíbrio seja naturalmente atingida. O desenvolvimento das fronteiras reativas acontece, naturalmente, quando o SPH co- meça a ser aplicado a problemas em superfície livre (MONAGHAN, 1994). Inicialmente, a força da parede assume a forma: �fik = K (( Δx ‖�rik‖ )p1 − ( r0 ‖�rik‖ )p2 ) �rik ‖�rik‖2 (64) sendo que p1 e p2 devem satisfazer a condição p1 > p2. Para os valores típicos, são reportados (MONAGHAN, 1994) p1 = 4 e p2 = 2 e K é escolhido conforme o tipo de problema. Para o caso de ruptura de barragens, toma-se K = gH0. Posteriormente, Monaghan e Kos (1999) descrevem a força exercida pela parede como uma composição de duas parcelas, uma vertical (Q) e uma horizontal (P), sendo: �fik = P(ψik)Q(ξik)�n (65) 54 Capítulo 4. O MÉTODO SPH A figura 6 ilustra uma situação com o emprego das fronteiras reativas. Figura 6 – Croqui esquemático com as distâncias horizontais (ψ) e verticais (ξ ) entre as partí- culas da fronteira e fluida. Fonte: Do próprio autor. Já Monaghan, Kos e Issa (2003) alteram os valores dos componentes Q e P propostos em Monaghan e Kos (1999), inserindo uma formulação que varia com a distância, semelhante ao comportamento do núcleo de suavização. Uma nota importante e inovadora das fronteira re- ativas apresentadas por Monaghan, Kos e Issa (2003) é o tratamento do escoamento nos cantos da fronteira. Os autores modificam o vetor normal se a partícula fluida está na região entre as normais das paredes que compõe o canto. O vetor normal, nesse caso, será dado pelo vetor uni- tário que liga o canto à partícula fluida. Notam-se, pela alteração do vetor normal e a utilização de uma parcela na composição da força da parede semelhante a um núcleo de suavização, os preconizadores da formulação de fronteira radial, apresentada em Kajtar e Monaghan (2008). Finalmente, Kajtar e Monaghan (2008) alteram a forma de como a força de fronteira é calculada. Inicialmente, os autores mostram, através da transformada de Fourier da força da fronteira, que esta é proporcional à função do núcleo de suavização. Posteriormente, os autores mostram que a força tangencial tradicionalmente inserida nos modelos têm baixa intensidade se comparada com a força radial. Sendo assim, Kajtar e Monaghan (2008) desprezam os termos tangenciais da força de fronteira (P na equação 65) e chegam a: �fik = 0,01 βm c2 i�rik W (rik/h) (66) 4.11 INSTABILIDADES NO SPH De modo geral, podem surgir algumas instabilidades no SPH, um delas sendo caracte- rizada por um aglutinamento de partículas que não tem correspondência com nenhum aspecto físico. Esse aglutinamento acontece devido a pressões negativas, induzidas pela equação de es- tado, por isso é conhecida como instabilidade de tensão (tensile instability). Um dos primeiros 4.12. EVOLUÇÃO TEMPORAL 55 trabalhos que identificaram a instabilidade de tensão foi Swegle, Hicks e Attaway (1995), atra- vés de uma análise para definir os critérios de estabilidade do SPH. Os autores utilizaram como referências o estado de tensão no escoamento e a derivada segunda do núcleo de suavização. Apesar de apontar as fontes da instabilidade, que são basicamente o aparecimento de pressões negativas, os autores não propõe um método para a sua correção, além de recomendar cuidados adicionais na escolha do núcleo de suavização. Seguiram-se posteriormente estudos a respeito da estabilidade do SPH, tal como Morris (1996). Mais tarde, Monaghan (2000) propõe que a instabilidade de tensão pode ser eliminada a partir da inserção de uma tensão artificial. A ideia é garantir que o aglutinamento seja evitado, a partir de um termo adicional de repulsão entre partículas. Esse termo seria então somado ao balanço de quantidade de movimento (equação 60), e representa a repulsão física entre os átomos em um sólido real. No caso de escoamentos de fluidos, essa tensão artificial se reduz a uma pressão artificial, dada por: Ri j = ( e1 |pi| ρi ρ j + e2 |p j| ρi ρ j )( W (ri j) W (Δx) )4 , (67) sendo e1 = 0,01 caso pi > 0, caso contrário e1 = 0,2. Da mesma forma, e2 = 0,01 caso p j > 0, caso contrário e2 = 0,2. A equação de quantidade de movimento fica então, sem considerar as condições de fronteira e as forças de corpo: d�ui dt =−∑ j m j ( pi + p j ρi ρ j +Πi j +Ri j ) �∇iWi j (68) 4.12 EVOLUÇÃO TEMPORAL Considerando, inicialmente, que o sistema de equações do SPH pode ser escrito na forma: d�r dt = �u (69) d�u dt = �F (70) dρ dt = D (71) e no caso de considerar a variação do comprimento de suavização h: dh dt = H . (72) Para solucionar esse sistema de equações explícitas, pode ser aplicado qualquer mé- todo de integração de equações diferenciais de primeira ordem, como o Leap-Frog ou Preditor- Corretor, desde que a máxima ordem do erro do método seja O(Δh2), sendo Δh o incremento 56 Capítulo 4. O MÉTODO SPH do método. A limitação com relação à ordem do erro se dá, principalmente, pelo erro intrínseco do SPH. Sabe-se, por expansão em séries de Taylor, que o SPH possui um erro da ordem de O(h2), portanto não há razão para utilizar métodos de integração de ordem superior a O(Δh2). Monaghan (2005) cita que as propriedades conservativas do método de integração são mais atrativas do que uma elevada ordem de precisão. Para apoiar essa assertiva, Monaghan (2005) fez uma simulação utilizando dois métodos de integração: o Runge-Kutta de 4a ordem e um integrador do tipo Verlet. Embora seja de segunda ordem, os resultados obtidos com o integrador Verlet são melhores que aqueles obtidos com o esquema Runge-Kutta. A justifica- tiva reside no fato de que o integrador Verlet conserva a quantidade de movimento angular e é reversível. Faz-se importante ressaltar, entretanto, que há certa relação entre as condições de contorno e os métodos de integração. Nas simulações feitas por Monaghan (2005), foram utilizadas as condições de contorno do tipo fronteira reativa. Ilustram-se, brevemente, os passos de alguns dos métodos de integração utilizados 4.12.1 Preditor-Corretor O método Preditor-Corretor, ou Euler Melhorado, pode ser sintetizado nos seguintes passos: a) Calcular�u, �F e D no instante atual (t); b) Calcular�u, �F e D no instante intermediário (t +Δt/2); c) Os valores em t +Δt serão dados pela equação do método: �r t+Δt = �r t + 1 2 Δt [�ut +�ut+Δt/2] (73) �u t+Δt = �u t + 1 2 Δt [ �F t + �F t+Δt/2] (74) ρ t+Δt = ρ t + 1 2 Δt [D t +D t+Δt/2] (75) 4.12.2 Método Verlet O método de integração simplético Verlet utiliza-se das características conservativas do integrador simplético (KAJTAR; MONAGHAN, 2008). Este método pode ser resumido nos seguintes passos: a) Calcular �F t e D t (t é o instante atual); b) Calcular �r t+Δt/2 e ρ t+Δt/2 no instante intermediário, fazendo: �r t+Δt/2 = �r t + 1 2 Δt�u t , (76) ρ t+Δt/2 = ρ t + 1 2 ΔtD t ; (77) 4.12. EVOLUÇÃO TEMPORAL 57 c) Com as coordenadas no tempo intermediário definidas, pode-se proceder o cálculo de �F t+Δt/2; d) A última etapa consiste em determinar �u t+Δt e �r t+Δt , fazendo: �u t+Δt = �u t +Δt �F t+Δt/2, (78) �r t+Δt = �r t+Δt/2 + 1 2 Δt�u t+Δt ; (79) e) Com �u t+Δt e �r t+Δt , calcula-se D t+Δt ; f) Finalmente, determina-se ρ t+Δ fazendo: ρ t+Δ = ρ t+Δt/2 + 1 2 ΔtD t+Δt . (80) 4.12.3 Runge-Kutta O método de Runge-Kutta de segunda ordem é empregado em algumas simulações, ao longo deste trabalho. O método pode ser descrito como: a) Calcular �F t e D t no instante atual; b) Calcular�rt+Δt/2,�ut+Δt/2 e ρ t+Δt/2 no instante intermediário; �r t+Δt/2 = �r t + 1 2 Δt�ut , (81) �u t+Δt/2 = �u t + 1 2 Δt �F t , (82) ρ t+Δt/2 = ρ t + 1 2 Δt D t ; (83) c) Com os valores intermediários, calcular �F t+Δt/2 e D t+Δt/2: d) Finalmente, calcular as grandezas no tempo final: �r t+Δt = �r t + 1 2 Δt�ut+Δt/2, (84) �u t+Δt = �u t + 1 2 Δt �F t+Δt/2, (85) ρ t+Δt = ρ t + 1 2 Δt D t+Δt/2; (86) As equações de evolução temporal valem para todas as partículas fluidas. Para casos especiais, por exemplo, na presença de um sólido com movimento livre, devem ser observadas restrições adicionais, que serão tratadas oportunamente (seção 7.2.3, p. 95). As condições de estabilidade numérica também variam de problema para problema. 58 Capítulo 4. O MÉTODO SPH 4.13 COMPARAÇÃO ENTRE MÉTODOS NUMÉRICOS: SPH e MPS É interessante notar que, embora seja um método que reproduza bem a incompressibi- lidade do escoamento, a função núcleo do MPS não possui derivada primeira contínua, o que resulta em maior oscilação para um mesmo número de partículas em comparação com o SPH (MONAGHAN; KOS; ISSA, 2003). Há, também, uma variante do SPH que descarta a necessidade de uma equação de estado para a pressão, assim como no MPS. Shao (2006) usa essa variante do SPH para analisar a quebra da onda e galgamento em barragens, incluindo efeitos turbulentos em sua modelagem. Diante do exposto, decidiu-se utilizar o SPH como método numérico neste trabalho. 59 5 ESCOAMENTO DE FLUIDOS IDEAIS De modo a validar o código apresentado na Parte I, apresentam-se testes clássicos da literatura para geração de ondas e escoamentos em superfície livre. Nesta fase, serão avalia- dos os casos de ruptura de barragens (dam break), a geração de ondas por meio de batedor (wavemaker) e o esvaziamento de um reservatório, através de um orifício. Nos problemas ali- nhavados, observa-se principalmente o comportamento do código no que tange à representação de escoamentos ideais bidimensionais. 5.1 RUPTURA DE BARRAGENS O problema típico de ruptura de barragem ou barreira é caracterizado por uma porção de líquido retido por uma comporta, que subitamente é removida. A massa liquida avança, en- tão, como uma onda, com velocidade da frente definida. Tal categoria de fenômeno apresenta grande interesse do ponto de vista de Engenharia, e devido à facilidade de implementação, é comumente usado como um caso teste para verificação de códigos numéricos em desenvolvi- mento. A figura 7 ilustra um ensaio experimental do problema em questão. Figura 7 – À esquerda, um croqui esquemático de uma ruptura de barragens, enquanto que à direita, observa-se o líquido retido, em ensaio experimental. Fonte: Do próprio autor. Fonte: Minussi (2007). Utilizam-se, neste primeiro teste, os resultados obtidos por Martin e Moyce (1952). Os autores usaram um canal de pequenas dimensões, 5,72×5,72×12,7 cm, de modo que a altura e a largura do fluido retido pudessem ser alteradas, de modo a possibilitar a obtenção de diversas relações geométricas. 60 Capítulo 5. ESCOAMENTO DE FLUIDOS IDEAIS 5.1.1 Descrição do modelo numérico Representam-se as variáveis do problema de ruptura de barragens através da altura do material retido, tal como ilustrado na figura 7. O material retido é modelado como um fluido ideal, com massa específica ρ = 1000 kg/m3. A lâmina d’água varia nos ensaios em dois valo- res: 2,86 e 5,72 cm. Os efeitos da viscosidade artificial estão presentes, com α = 0,3 e β = 0. O núcleo de suavização usado é o spline cúbico 2D, cujo comprimento de suavização vale h = 1,2Δx, sendo que Δx variou entre 1,25 e 2,5 mm, com o propósito de manter a resolu- ção constante (ou seja, o número de partículas distribuídas na vertical). A celeridade do som da simulação vale c = 23 m/s. As condições de contorno são representadas pelas fronteiras reativas, sendo que a força na parede é dada pelo produto de um termo horizontal (P) e outro vertical (Q, apresentados na equação 65), são dados por: P(ψ) = { 1− ψ Δp , se 0 < ψ < Δp 0, se 0 ≥ ψ ≥ Δp (87) Q(ξ ) = ⎧⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎩ 2 3 K, se 0 < ξ/h < 2/3 K ( 2 ξ h − 3 2 ( ξ h )2 ) , se 2/3 < ξ/h < 1 1 2 K (2−ξ/h)2, se 1 < ξ/h < 2 0, se ξ ≥ 2 (88) sendo K = 0,02c2/ξ . Como condição inicial, foi imposto o gradiente de pressão teórico para todas as partículas fluidas. Usa-se a correção XSPH para deslocamento das partículas, com ε = 0,5, e a clássica equação de estado para a pressão (equação 56). O passo de tempo Δt obedece às restrições impostas pela condição CFL (Courant- Friedrichs-Lewy) e é o menor valor entre Δt1 e Δt2, sendo Δt1 = hi/(c+ 0,6α) calculado para todas as partículas fluidas e Δt2 = y/(0,1c), calculado para todas as partículas interagindo com a fronteira. O integrador numérico utilizado foi do tipo Preditor-Corretor (apresentado na seção 4.12.1). 5.1.2 Resultados De modo a facilitar na organização, os resultados são agrupados e apresentados com mesma relação H0/L0. Sendo assim, para o caso de H0/L0 = 1, foram escolhidas duas configu- rações: uma com H0 = 2,86 cm e outra com H0 = 5,72 cm. No caso de H0 = 2,86 cm, tem-se a respectiva simulação SPH, contando com quase 1500 partículas para representação do fluido escoante e fronteiras do problema. A figura 8a mostra a comparação dos resultados obtidos com os resultados experimentais de Martin e Moyce (1952). Já na figura 8b, faz-se a mesma compa- ração, mas para o segundo caso, com H0 = 5,72 cm, contando com uma simulação com pouco mais de 1500 partículas. 5.1. RUPTURA DE BARRAGENS 61 Figura 8 – Comparação numérico-experimental, para uma relação H0/L0 = 1 com (a) H0 = 2,86 cm e (b) H0 = 5,72 cm. Fonte: Do próprio autor. Para o caso H0/L0 = 2, foram também ensaiados para duas composições geométricas: H0 = 5,72 cm e H0 = 2,86 cm. No primeiro caso, foram geradas partículas com espaçamento médio de 1,25 mm, o que resultou em aproximadamente de 2800 partículas. Mantendo o mesmo Δx, para a segunda altura de material retido, obteve-se pouco mais de 1400 partículas. Os resul- tados, com a respectiva simulação numérica, são mostrados nas figuras 9a e 9b. Figura 9 – Comparação numérico-experimental, para uma relação H0/L0 = 2 com (a) H0 = 5,72 cm e (b) H0 = 2,86 cm. Fonte: Do próprio autor. Nas figuras 8 e 9, os resultados numéricos foram normalizados, sendo Z = z/L0 o al- cance adimensional da frente, e T = t √ H0 g/L2 0 o tempo adimensional. Tanto nos dados ex- perimentais quanto nos numéricos, há um ajuste no zero, mediante a incapacidade do aparato experimental na determinação precisa do tempo inicial do ensaio. Em termos práticos, esse procedimento apenas translada a envoltória obtida para o alcance da frente. 62 Capítulo 5. ESCOAMENTO DE FLUIDOS IDEAIS É possível notar em algumas figuras uma repentina mudança de declividade na ten- dência dos resultados experimentais. Esse fato pode ser explicado pela influência dos aspectos viscosos do escoamento, uma vez que as dimensões características do ensaio e notadamente da frente da ruptura são relativamente pequenas. No entanto, a análise de fluido ideal restringe-se à etapa inercial da ruptura, evidenciada nos primeiros momentos do fenômeno. Nesta região, principalmente nos casos em que a frente de ruptura possui menor razão H0/L0 (figuras 8a e 8b), há concordância entre resultados experimentais e numéricos. Essas verificações comprovam a validade do modelo para ruptura de barragem em fluido ideal. 5.2 GERAÇÃO, PROPAGAÇÃO E QUEBRA DE ONDAS É sabido que o padrão de ondas incidente em áreas costeiras é de fundamental impor- tância, não só pela determinação de zonas com baixa influência das ondas, mas também para o dimensionamento de obras de contenção e proteção. A dinâmica das ondas em áreas costeiras, objeto de estudo da Hidráulica Marítima, está assentada, de modo geral, na existência de um potencial de velocidades donde derivam as propriedades cinemáticas do escoamento. Uma aproximação analítica apresenta dificuldades, uma vez que as equações que regem o problema da onda, com diminuição da profundidade, têm forte característica não-linear. Nesse contexto, avaliar ainda a dinâmica da quebra da onda, torna o problema ainda mais complexo