Latest revision as of 09:58, 10 January 2018

Resumo

Este artigo apresenta um novo esquema explícito para a solução das equações de águas rasas em uma e duas dimensões, desenvolvido a partir do método dos elementos de conservação e elementos de solução espaço-tempo, aqui abreviado por método CE/SE. As funções de base utilizadas são expansões de Taylor de segunda ordem no tempo e no espaço. Esse aumento na ordem das funções de aproximação produz o aumento no número de variáveis de marcha no tempo, por isso, além das variáveis fluxo e suas inclinações, também são incógnitas no presente esquema suas derivadas espaciais de segunda ordem. Um processo iterativo para o cálculo das derivadas de primeira e segunda ordem é formulado para problemas com choques e descontinuidades. Experimentos computacionais demonstram acurácia de terceira ordem. Os problemas de ruptura de barragem unidimensional e bidimensional considerados validam a acurácia e robustez do esquema.

Palavras-Chave: Leis de conservação, Equações de águas rasas, Volume de controle espaço-tempo, Alta ordem de acurácia.

Abstract

This paper presents a new explicit scheme for the solution of shallow water equations in one and two space dimensions, developed from the space-time conservation element and solution element (CE/SE) method. The basis functions used are second-order Taylor expansions in time and space. This increase in the order of the approximation functions produces an increase in the number of unknowns in the scheme, therefore, besides the flow variables and their slopes, their second-order partial derivatives are also unknown in the present scheme. An iterative process for the calculation of the first and second order derivatives is formulated for problems with shocks and discontinuities. Computational experiments demonstrate third-order accuracy. The one-dimensional and two-dimensional dam-break problems presented validate the accuracy and robustness of this scheme.

1 Introdução

As equações de águas rasas constituem um dos modelos mais comumente usados na análise de fluxo de água em rios ou áreas costais [1], sobretudo no que tange à simulações de fluxo de ruptura de barragem [2,3,4,5], parte essencial do projeto e avaliação da segurança de barragens, controle de inundação de rios, mitigação de desastres de bacias hidrográficas, etc [6]. Além da importância prática em hidráulica e engenharia costal, fornecem um ótimo modelo matemático para equações diferenciais hiperbólicas não lineares que podem ter soluções como ondas de choque [7] e, por isso, também são frequentemente utilizadas como problema de referência para novos esquemas de previsão numérica [8,9].

Extensivos estudos numéricos foram feitos no sentido de analisar os fenômenos governados por essas equações e diversos métodos foram desenvolvidos, como de volumes finitos [10,11,12], elementos finitos [13,14,15,16], diferenças finitas [17,18] e outros [19,20].

Recentemente, um novo método baseado em volumes de controle, denominado método dos elementos de conservação e elementos de solução espaço-tempo, do inglês, space-time conservation element and solution element method (CE/SE), foi proposto por Chang e To (1991) [21] e aplicações do mesmo em problemas de ruptura de barragem logo foram feitas [22,23,24]. Este método foi desenvolvido com o objetivo de obter soluções numéricas de leis de conservação e possui características importantes [25,26,27,28]: sua construção combina informações de ambas as formas diferencial e integral das leis de conservação; o conjunto de variáveis de marcha é formado pelas variáveis fluxo e suas derivadas espaciais; não são utilizadas técnicas de interpolação ou extrapolação nos valores da malha; não são utilizadas técnicas baseadas nas características, o que torna o método de fácil extensão à várias dimensões, bem como aplicável à equações não hiperbólicas, como as equações de Navier-Stokes, por exemplo; possui uma simples notação estêncil, o que facilita a programação. É caracterizado, também, por sua versatilidade, tendo em vista que pode ser adaptado à malhas não uniformes e não estruturadas [29,30], bem como pode ser estendido à altas ordens [31,32], razão pela qual possui aplicações nas mais diversas áreas [33,34,35,36].

O objetivo deste trabalho é apresentar a construção de um novo esquema CE/SE explícito de alta ordem para as equações de águas rasas e verificar sua eficiência em problemas de ruptura de barragem. Para isso, este artigo está organizado da seguinte maneira: na seção 2 é apresentado o modelo matemático a ser estudado (subseção 2.1), bem como é desenvolvido o esquema numérico discreto (subseção 2.2); na seção 3, experimentos numéricos unidimensionais e bidimensionais são considerados com a finalidade de validar o método; as conclusões são tecidas junto à seção 4 e, ao final, na seção 5, dispõe-se a lista de referências.

2 Equações Governantes e Método CE/SE

2.1 Equações de Águas Rasas

As equações de águas rasas, também conhecidas como equações de Saint Venant, possuem a seguinte forma vetorial conservativa

em que

sendo ${\textstyle h}$ a profundidade do fluxo (m); ${\textstyle u}$ e ${\textstyle v}$ são as velocidades médias do escoamento (m/s) nas direções ${\textstyle x}$ e ${\textstyle y}$, respectivamente; ${\textstyle g}$ é a aceleração da gravidade (${\textstyle 9.812}$m/s${\textstyle ^{2}}$); ${\textstyle S_{0{x}}=\partial Z_{b}/\partial x}$ e ${\textstyle S_{0{y}}=\partial Z_{b}/\partial y}$ são as inclinações do fundo do canal cuja topografia é dada por ${\textstyle Z_{b}}$; ${\textstyle S_{f{x}}}$ e ${\textstyle S_{f{y}}}$ são os termos de resistência ao escoamento, nas direções ${\textstyle x}$ e ${\textstyle y}$, respectivamente, definidos por

em que ${\textstyle n}$ (s/m${\textstyle ^{1/3}}$) é o coeficiente de fricção de Manning.

2.2 Esquema CE/SE

Considere, por simplicidade, a Eq. (1) como

então, pelo teorema da divergência no espaço ${\textstyle \mathbb {R} ^{3}}$, tem-se que a Eq. (5) representa a forma diferencial da lei integral de conservação

em que ${\textstyle {\textbf {H}}_{m}=(f_{m},g_{m},q_{m})}$, ${\textstyle m=1,2,3}$ e ${\textstyle S(V)}$ representa o contorno de uma região espaço-tempo ${\textstyle V\subset \mathbb {R} ^{3}}$.

No método CE/SE, um elemento de conservação (CE) é uma região espaço-tempo em que a conservação do fluxo, Eq. (6), é forçada, enquanto que um elemento de solução (SE) é uma região espaço-tempo, normalmente distinta, em que as variáveis fluxo são supostamente suaves e a Eq. (5) é válida. Para definir estes dois importante objetos, considere primeiramente uma malha no plano ${\textstyle x}$-${\textstyle y}$, conforme Fig. 1a. Existem dois grupos de pontos, marcados por círculos e cruzes, que representam nós da malha em dois níveis de tempo diferentes. A cada ponto ${\textstyle (i,j,k)}$ da malha associa-se um CE e um SE. O CE é definido como o quadrilátero ${\textstyle EFGHE'F'G'H'}$ e o SE é a união do quadrilátero ${\textstyle P''Q''R''S''P'Q'R'S'}$ e o polígono ${\textstyle EFGH}$ (veja Fig. 1b).

(a)	(b)
Figura 1: (a) Pontos da malha representativa no plano ${\displaystyle x}$-${\displaystyle y}$ e (b) as definições dos elementos de solução SE${\displaystyle (i,j,k)}$ e elementos de conservação CE${\displaystyle (i,j,k)}$ no ponto ${\displaystyle (i,j,k)}$ da malha [23,37].

Para todo ${\textstyle (x,y,t)\in {\mbox{SE}}(i,j,k)}$, aproxima-se ${\textstyle q_{m}}$ por um polinômio de Taylor de segunda ordem

Aproximações análogas são feitas sobre as funções ${\textstyle f_{m}}$ e ${\textstyle g_{m}}$, ${\textstyle m=1,2,3}$. O termo fonte, por outro lado, é aproximado por uma ordem a menos:

Dessa forma, a Eq. (5) pode ser aproximada no SE${\textstyle (i,j,k)}$ por div ${\textstyle {\textbf {H}}_{m}^{*}=S_{m}^{*}}$, em que ${\textstyle {\textbf {H}}_{m}^{*}=\left(f_{m}^{*},g_{m}^{*},q_{m}^{*}\right)}$, ${\textstyle m=1,2,3}$, ao mesmo tempo em que a Eq. (6) é aproximada no CE${\textstyle (i,j,k)}$ por

2.2.1 Avaliação de (q_m)_i,j^k

Substituindo as funções ${\textstyle q_{m}^{*},f_{m}^{*},g_{m}^{*}}$ e ${\textstyle S_{m}^{*}}$ em (9), após diversas simplificações, obtém-se o esquema de avanço no tempo para a variável ${\textstyle (q_{m})_{i,j}^{k}}$:

em que

É importante notar que as Equações (10) e (11) dependem apenas das incógnitas ${\textstyle q_{m},(q_{m})_{x},(q_{m})_{y},(q_{m})_{xx},(q_{m})_{yy}}$ e ${\textstyle (q_{m})_{xy}}$ do tempo ${\textstyle t_{k-1/2}}$, para ${\textstyle m=1,2,3}$ pois, pelas Equações (1)-(4), ${\textstyle f_{m}}$, ${\textstyle g_{m}}$ e ${\textstyle S_{m}}$ também dependem destas últimas. Por outro lado, é preciso conhecer ${\textstyle [(q_{m})_{xx}]_{i,j}^{k}}$ e ${\textstyle [(q_{m})_{yy}]_{i,j}^{k}}$ previamente à obtenção de ${\textstyle (q_{m})_{i,j}^{k}}$ no tempo ${\textstyle t_{k}}$, conforme Eq. (10). A seção a seguir descreve a obtenção dessas derivadas espaciais duplas.

2.2.2 Avaliação de [(q_m)_xx]_i,j^k, [(q_m)_xy]_i,j^k e [(q_m)_yy]_i,j^k

Para prosseguir com o cálculo de ${\textstyle [(q_{m})_{xx}]_{i,j}^{k}}$, constrói-se primeiramente uma equação auxiliar a partir da Eq. (5), derivando-a duas vezes em relação a ${\textstyle x}$, obtendo as seguintes formas diferencial e integral

As Equações em (12) possuem, conforme Equações (7) e (8), para todo ${\textstyle (x,y,t)\in }$ SE${\textstyle (i,j,k)}$, os seguintes análogos numéricos:

O campo vetorial aproximado ${\textstyle {\hat {\textbf {H}}}_{m}^{*}=((f_{m}^{*})_{xx}-(S_{m}^{*})_{x},(g_{m}^{*})_{xx},(q_{m}^{*})_{xx})}$ é constante, de modo que a avaliação da Eq. (13) sobre o elemento de conservação descrito pela Fig. 1a retorna, para ${\textstyle m=1,2,3}$:

Para computar as variáveis de marcha ${\textstyle [(q_{m})_{yy}]_{i,j}^{k}}$ e ${\textstyle [(q_{m})_{xy}]_{i,j}^{k}}$ procede-se de modo análogo, obtendo-se após todos os cálculos e simplificações:

Aborda-se na seção a seguir o cálculo das derivadas de primeira ordem das variáveis dinâmicas ${\textstyle q_{m}}$, ${\textstyle m=1,2,3}$.

2.2.3 Avaliação de [(q_m)_x]_i,j^k e [(q_m)_y]_i,j^k

As variáveis ${\textstyle [(q_{m})_{x}]_{i,j}^{k}}$ e ${\textstyle [(q_{m})_{y}]_{i,j}^{k}}$ podem ser determinadas por meio de uma estratégia análoga àquela apresentada na seção [[#2.2.2 Avaliação de ${\displaystyle [(q_{m})_{xx}]_{i,j}^{k}}$, ${\displaystyle [(q_{m})_{xy}]_{i,j}^{k}}$ e ${\displaystyle [(q_{m})_{yy}]_{i,j}^{k}}$|2.2.2]] anterior, bastando, para isso, construir leis diferenciais de conservação (já nas formas aproximadas)

e integrais

As avaliações das equações em (18) fornecem as inclinações equacionadas por

para ${\textstyle \xi =x,y}$ e ${\textstyle m=1,2,3}$, em que

O esquema formado pelo conjunto de Equações (10), (14)-(16) e (19) é interessante apenas em problemas com soluções de comportamento suave. Para problemas com descontinuidades ou formação de choque, é necessário modificar o cálculo das variáveis em (14)-(16) e (19). A proposta apresentada neste artigo é uma variação do procedimento adotado por [38,39] e consiste num processo iterativo.

Para prosseguir, considere as seguintes formas regressivas e progressivas:

tal que o operador ${\textstyle \delta _{\xi }^{\pm }}$ seja definido como

As aproximações para as derivadas serão definidas como

e a função ${\textstyle W}$ [26,24]:

Observe que neste caso, as Equações (24) e (25) representam médias ponderadas das diferenças progressivas e regressivas descritas em (22) e (23). No caso de ${\textstyle \alpha =0}$, (24) e (25) tornam-se diferenças finitas centrais [23]. Adota-se neste trabalho ${\textstyle \alpha =1}$.

Importante notar, no entanto, que as Equações (24) e (25) dependem da variável ${\textstyle (q_{m})_{i,j}^{k}}$ que, por sua vez, depende de ${\textstyle [(q_{m})_{xx}]_{i,j}^{k}}$ e ${\textstyle [(q_{m})_{yy}]_{i,j}^{k}}$, conforme Eq. (10). Dito de outra forma, para obter ${\textstyle (q_{m})_{i,j}^{k}}$ é necessário obter antes ${\textstyle [(q_{m})_{xx}]_{i,j}^{k}}$ e ${\textstyle [(q_{m})_{yy}]_{i,j}^{k}}$, e vice-versa. Para sanar esta dificuldade, construímos o processo iterativo descrito no Algoritmo 1.

A estabilidade do esquema CE/SE satisfaz a condição de Courant-Friedrichs-Lewy (CFL) que, para as equações de Saint Venant unidimensional e bidimensional são as restrições [23]

respectivamente, onde o número de Courant satisfaz ${\textstyle 0<CFL\leq 1}$. Maiores detalhes sobre estabilidade do método CE/SE pode ser encontrado em [40].

3 Exemplos Numéricos

Esta seção tem por objetivo avaliar e validar o esquema numérico desenvolvido. Para isso, alguns problemas teste são avaliados, como testes de acurácia em uma e duas dimensões, bem como problemas clássicos de ruptura de barragem unidimensionais (subseção 3.1) e bidimensionais (subseção 3.2).

3.1 Exemplos Unidimensionais

Exemplo 3.1.1: O objetivo deste exemplo teste é verificar experimentalmente a ordem de acurácia do esquema CE/SE. Considere o sistema hiperbólico linear unidimensional

e ${\displaystyle c}$ é uma constante dada. Impondo as condições inicias ${\displaystyle \phi (x,0)=-c^{-1}\sin(2\pi x)}$ e ${\displaystyle u(x,0)=0}$, é possível obter a seguinte solução exata [41,42]:

As soluções são computadas com as Equações (10), (14)-(16) e (19). A Tabela 1 apresenta os erros calculados nas normas ${\displaystyle L^{1}}$ e ${\displaystyle L^{2}}$, nas simulações realizadas sobre o domínio computacional ${\displaystyle -1\leq x\leq 1}$, ${\displaystyle c=1}$, ${\displaystyle CFL=0.8}$ até atingir o tempo ${\displaystyle t=1.2}$s, sendo ${\displaystyle n}$ o número de células no espaço. O experimento verifica a terceira ordem de acurácia do esquema para ambas as variáveis ${\displaystyle u}$ e ${\displaystyle \phi }$.