High-order CE/SE scheme for dam-break flow simulation

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

(1)

em que

{\textbf {Q}}=\left[{\begin{array}{c}q_{1}\\q_{2}\\q_{3}\end{array}}\right]=\left[{\begin{array}{c}h\\hu\\hv\end{array}}\right],\quad {\textbf {F}}=\left[{\begin{array}{c}f_{1}\\f_{2}\\f_{3}\end{array}}\right]=\left[{\begin{array}{c}hu\\hu^{2}+gh^{2}/2\\huv\end{array}}\right]=\left[{\begin{array}{c}q_{2}\\q_{2}^{2}/q_{1}+gq_{1}^{2}/2\\q_{2}q_{3}/q_{1}\end{array}}\right],

(2)

{\textbf {G}}=\left[{\begin{array}{c}g_{1}\\g_{2}\\g_{3}\end{array}}\right]=\left[{\begin{array}{c}hv\\huv\\hv^{2}+gh^{2}/2\end{array}}\right]=\left[{\begin{array}{c}q_{3}\\q_{2}q_{3}/q_{1}\\q_{3}^{2}/q_{1}+gq_{1}^{2}/2\end{array}}\right]\quad {\textrm {e}}\quad {\textbf {S}}=\left[{\begin{array}{c}S_{1}\\S_{2}\\S_{3}\end{array}}\right]=\left[{\begin{array}{c}0\\gq_{1}(S_{0x}-S_{fx})\\gq_{1}(S_{0y}-S_{fy})\end{array}}\right],

(3)

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

S_{fx}=\displaystyle {\frac {n^{2}u{\sqrt {u^{2}+v^{2}}}}{h^{4/3}}}={\frac {n^{2}q_{2}{\sqrt {q_{2}^{2}+q_{3}^{2}}}}{q_{1}^{10/3}}}\quad {\textrm {e}}\quad S_{fy}=\displaystyle {\frac {n^{2}v{\sqrt {u^{2}+v^{2}}}}{h^{4/3}}}={\frac {n^{2}q_{3}{\sqrt {q_{2}^{2}+q_{3}^{2}}}}{q_{1}^{10/3}}},

(4)

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

(5)

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

(6)

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 $x$ - $y$ e (b) as definições dos elementos de solução SE $(i,j,k)$ e elementos de conservação CE $(i,j,k)$ no ponto $(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

(7)

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:

(8)

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

(9)

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}}$ :

(q_{m})_{i,j}^{k}+\displaystyle {\frac {\Delta x^{2}}{48}}[(q_{m})_{xx}]_{i,j}^{k}+{\frac {\Delta y^{2}}{48}}[(q_{m})_{yy}]_{i,j}^{k}=[P_{m}^{y}(\Delta x,f_{m},g_{m})]_{i-1/2,j}^{k-1/2}+[P_{m}^{y}(-\Delta x,f_{m},g_{m})]_{i+1/2,j}^{k-1/2}+

(10)

em que

(11)

É 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

\displaystyle {\frac {\partial (q_{m})_{xx}}{\partial {t}}}+{\frac {\partial }{\partial x}}\left[(f_{m})_{xx}-(S_{m})_{x}\right]+{\frac {\partial (g_{m})_{xx}}{\partial y}}={\textrm {div}}{\hat {\textbf {H}}}_{m}=0\quad {\textrm {e}}\quad \oint _{S(V)}{\hat {\textbf {H}}}_{m}\cdot d{\textbf {s}}=0,\quad m=1,2,3.

(12)

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:

(13)

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}$ :

(14)

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:

${[}(q_{m})_{yy}{]}_{i,j}^{k}={\frac {1}{4}}\left\{[(q_{m})_{yy}]_{i-1/2,j}^{k-1/2}+[(q_{m})_{yy}]_{i+1/2,j}^{k-1/2}+[(q_{m})_{yy}]_{i,j-1/2}^{k-1/2}+\right.$ $\left.[(q_{m})_{yy}]_{i,j+1/2}^{k-1/2}\right\}+{\frac {\Delta t}{2\Delta y}}\left\{[(g_{m})_{yy}]_{i,j-1/2}^{k-1/2}\right.$
$\left.-[(g_{m})_{yy}]_{i,j+1/2}^{k-1/2}+[(S_{m})_{y}]_{i,j+1/2}^{k-1/2}-[(S_{m})_{y}]_{i,j-1/2}^{k-1/2}\right\}+$ ${\frac {\Delta t}{2\Delta x}}\left\{[(f_{m})_{yy}]_{i-1/2,j}^{k-1/2}-[(f_{m})_{yy}]_{i+1/2,j}^{k-1/2}\right\},$	(15)
${[}(q_{m})_{xy}{]}_{i,j}^{k}={\frac {1}{4}}\left\{[(q_{m})_{xy}]_{i-1/2,j}^{k-1/2}+[(q_{m})_{xy}]_{i+1/2,j}^{k-1/2}+[(q_{m})_{xy}]_{i,j-1/2}^{k-1/2}+[(q_{m})_{xy}]_{i,j+1/2}^{k-1/2}\right\}$
$+{\frac {\Delta t}{2\Delta x}}\left\{[(f_{m})_{xy}]_{i-1/2,j}^{k-1/2}-[(f_{m})_{xy}]_{i+1/2,j}^{k-1/2}+{\frac {1}{2}}[(S_{m})_{y}]_{i+1/2,j}^{k-1/2}-{\frac {1}{2}}[(S_{m})_{y}]_{i-1/2,j}^{k-1/2}\right\}$
$+{\frac {\Delta t}{2\Delta y}}\left\{[(g_{m})_{xy}]_{i,j-1/2}^{k-1/2}-[(g_{m})_{xy}]_{i,j+1/2}^{k-1/2}+{\frac {1}{2}}[(S_{m})_{x}]_{i,j+1/2}^{k-1/2}-{\frac {1}{2}}[(S_{m})_{x}]_{i,j-1/2}^{k-1/2}\right\},\quad m=1,2,3.$	(16)

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 $[(q_{m})_{xx}]_{i,j}^{k}$ , $[(q_{m})_{xy}]_{i,j}^{k}$ e $[(q_{m})_{yy}]_{i,j}^{k}$ |2.2.2]] anterior, bastando, para isso, construir leis diferenciais de conservação (já nas formas aproximadas)

(17)

e integrais

(18)

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

(19)

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

$Q_{m}^{(1)}(\Delta x,\xi )=\displaystyle {\frac {1}{4}}\left\{(q_{m})_{\xi }+{\frac {\Delta x}{4}}(q_{m})_{x\xi }+{\frac {\Delta t}{4}}\left[(f_{m})_{x\xi }-(g_{m})_{y\xi }\right]+\right.$ $\left.{\frac {2\Delta t}{\Delta x}}\left[(f_{m})_{\xi }+\displaystyle {\frac {\Delta t}{4}}(f_{m})_{t\xi }\right]+{\frac {\Delta t}{2}}(S_{m})_{\xi }\right\},$	(20)
$Q_{m}^{(2)}(\Delta y,\xi )=\displaystyle {\frac {1}{4}}\left\{(q_{m})_{\xi }+{\frac {\Delta y}{4}}(q_{m})_{y\xi }+{\frac {\Delta t}{4}}\left[(g_{m})_{y\xi }-(f_{m})_{x\xi }\right]+\right.$ $\left.{\frac {2\Delta t}{\Delta y}}\left[(g_{m})_{\xi }+\displaystyle {\frac {\Delta t}{4}}(g_{m})_{t\xi }\right]+{\frac {\Delta t}{2}}(S_{m})_{\xi }\right\}.$	(21)

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:

(22)

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

\delta _{\xi }^{\pm }(\phi )=\displaystyle \left\{\displaystyle {\begin{array}{cl}\displaystyle \pm {\frac {2}{\Delta \xi }}\left[(\phi ')_{i\pm {1}/2,j}^{k}-\phi _{i,j}^{k}\right],&{\textrm {se}}\;\xi =x\\\displaystyle \pm {\frac {2}{\Delta \xi }}\left[(\phi ')_{i,j\pm {1}/2}^{k}-\phi _{i,j}^{k}\right],&{\textrm {se}}\;\xi =y\end{array}}\right.\quad {\textrm {e}}\quad

\phi '=\left\{{\begin{array}{cl}\displaystyle \phi +{\frac {\Delta t}{2}}\phi _{t}+{\frac {\Delta t^{2}}{8}}\phi _{tt},&{\textrm {se}}\;\phi =(q_{m})_{x},(q_{m})_{y},\\\displaystyle \phi +{\frac {\Delta t}{2}}\phi _{t},&{\textrm {se}}\;\phi =(q_{m})_{xx},(q_{m})_{xy},(q_{m})_{yy}.\end{array}}\right.

(23)

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

$[(q_{m})_{\xi }]_{i,j}^{k}=W\left\{[(q_{m}^{+})_{\xi }]_{i,j}^{k},[(q_{m}^{-})_{\xi }]_{i,j}^{k},\alpha \right\},$	(24)
${[}(q_{m})_{\zeta \xi }{]}_{i,j}^{k}=W\left\{[(q_{m}^{+})_{\zeta \xi }]_{i,j}^{k},[(q_{m}^{-})_{\zeta \xi }]_{i,j}^{k},\alpha \right\},\quad$ $m=1,2,3,\quad \xi =x,y,\quad \zeta =x,y\quad {\textrm {e}}\quad \alpha \geq 0,$	(25)

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

W(x_{-},x_{+},\alpha )=\left\{{\begin{array}{cl}\displaystyle {\frac {|x_{+}|^{\alpha }x_{-}+|x_{-}|^{\alpha }x_{+}}{|x_{+}|^{\alpha }+|x_{-}|^{\alpha }}},&{\textrm {se}}\quad |x_{+}|^{\alpha }+|x_{-}|^{\alpha }\neq 0,\\0,&{\textrm {casocontrario.}}\end{array}}\right.

(26)

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.

Algoritmo 1: Algoritmo CE/SE para o cálculo das variáveis de marcha sobre um determinado ponto

da malha.

início

Defina uma tolerância

Faça

repita

Faça

;

Calcule

e

com as Equações (14) e (15), respectivamente;

Calcule

com a Eq. (10);

repita

Calcule

e

com a Eq. (24);

Calcule

,

e

com a Eq.(25);

Calcule

com a Eq. (10);

até

fim

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]

(27)

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

(28)

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

(29)

As soluções são computadas com as Equações (10), (14)-(16) e (19). A Tabela 1 apresenta os erros calculados nas normas $L^{1}$ e $L^{2}$ , nas simulações realizadas sobre o domínio computacional $-1\leq x\leq 1$ , $c=1$ , $CFL=0.8$ até atingir o tempo $t=1.2$ s, sendo $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 $u$ e $\phi$ .

Tabela 1. Ordem de acurácia experimental do método CE/SE para o sistema hiperbólico linear unidimensional.

	Norma $L^{1}$				Norma $L^{2}$
n	Erro ( $\phi$ )	Ordem	Erro ( $u$ )	Ordem	Erro ( $\phi$ )	Ordem	Erro ( $u$ )	Ordem
20	8.322 $\times 10^{-3}$	–	1.662 $\times 10^{-2}$	–	7.157 $\times 10^{-3}$	–	1.347 $\times 10^{-2}$	–
40	1.849 $\times 10^{-3}$	2.170	2.524 $\times 10^{-3}$	2.719	1.488 $\times 10^{-3}$	2.266	2.042 $\times 10^{-3}$	2.722
80	3.067 $\times 10^{-4}$	2.592	3.401 $\times 10^{-4}$	2.892	2.443 $\times 10^{-4}$	2.607	2.715 $\times 10^{-4}$	2.911
160	4.369 $\times 10^{-5}$	2.811	4.386 $\times 10^{-5}$	2.955	3.493 $\times 10^{-5}$	2.806	3.465 $\times 10^{-5}$	2.970
320	5.813 $\times 10^{-6}$	2.910	5.566 $\times 10^{-6}$	2.978	4.665 $\times 10^{-6}$	2.904	4.369 $\times 10^{-6}$	2.987
640	7.490 $\times 10^{-7}$	2.956	7.010 $\times 10^{-7}$	2.989	6.030 $\times 10^{-7}$	2.952	5.490 $\times 10^{-7}$	2.992

Exemplo 3.1.2: Considera-se agora um problema ideal de ruptura de barragem sobre um domínio molhado, isto é, a quebra de barragem é instantânea, o fundo é plano e não existe resistência ao escoamento. As condições iniciais para esta configuração seguem o clássico problema de Riemann

(30)

O domínio considerado é $0{\mbox{m}}$ $\leq x\leq 2000{\mbox{m}}$ , $x_{0}=1000{\mbox{m}}$ e a solução analítica pode ser encontrada em [43] ou [44]. Utiliza-se, para fins de análise de comportamento, a mesma relação de equações do Ex. 3.1.1 anterior, isto é, Equações (10), (14)-(16) e (19). As Fig. 2a e Fig. 2b demonstram o comportamento da solução numérica em relação à analítica, calculadas no tempo $t=52{\mbox{s}}$ , com $n+1=201$ pontos, incremento espacial $\Delta x=10{\mbox{m}}$ , $CFL=0.8$ e profundidades iniciais a montante e a jusante $h_{l}=10{\mbox{m}}$ e $h_{r}=5{\mbox{m}}$ , respectivamente. Note que a resposta numérica é coerente com a analítica, embora apresente suavidade que a distancie nas regiões com mudanças abruptas.


(a)	(b)
Figura 2. Elevação da superfície da água $h$ e velocidade $u$ para o problema de ruptura de barragem, computada no tempo $t=52{\mbox{s}}$ , utilizando uma malha uniforme com $201$ pontos, incremento espacial $\Delta x=10{\mbox{m}}$ e $h_{r}/h_{l}=0.5$ .

Utilizando o Algoritmo 1, sob as mesmas condições e com os mesmos parâmetros, obtém-se os gráficos constantes na Fig. 3. Observa-se que esta solução numérica é superior à anterior, sobretudo no que tange às regiões de rápidas mudanças. Conforme Zhang et al. (2012) [23], a razão $h_{r}/h_{l}$ é um importante índice para julgar a aplicabilidade e a acurácia de esquemas numéricos no modelo 1D de ruptura de barragem. Segundo os mesmos autores, os regimes de escoamentos subcrítico e supercrítico existem simultaneamente num canal sem fricção, horizontal e retangular, quando $h_{r}/h_{l}<0.138$ . Altera-se, neste sentido, estes parâmetros para uma razão $h_{r}/h_{l}=0.001$ . Os resultados simulados são mostrados nas Figs. 3c e 3d.


(a)	(b)

(c)	(d)
Figura 3. Elevação da superfície da água $h$ e velocidade $u$ para o problema de ruptura de barragem, no tempo $t=52{\mbox{s}}$ , numa malha uniforme com $201$ pontos e respectivo incremento espacial $\Delta x=10{\mbox{m}}$ , calculados com o Algoritmo 1, sendo: (a)-(b) $h_{l}=10{\mbox{m}}$ e $h_{r}=5{\mbox{m}}$ , (c)-(d) $h_{l}=10{\mbox{m}}$ e $h_{r}=0.01{\mbox{m}}$ .

3.2 Exemplos Bidimensionais

Exemplo 3.2.1: Considere agora o sistema hiperbólico linear bidimensional

(31)

e $c$ é uma constante dada. Para os dados iniciais $\phi (x,y,0)=-c^{-1}[\sin(2\pi x)+\sin(2\pi y)]$ e $u(x,y,0)=v(x,y,0)=0$ , este problema admite a seguinte solução exata [41,42]:

(32)

A Tabela 2 apresenta os erros numéricos calculados nas normas $L^{1}$ e $L^{2}$ , no tempo $t=0.2$ s, com parâmetros especificados em $c=1$ , $CFL=0.4$ e domínio computacional definido em $-1\leq x,y\leq 1$ , sendo $n_{x}\times n_{y}$ o número de células utilizadas na discretização. O experimento confirma novamente a acurácia de terceira ordem do esquema CE/SE.

Tabela 2. Ordem de acurácia experimental do método CE/SE para o sistema hiperbólico linear bidimensional..

	Norma $L^{1}$				Norma $L^{2}$
$n_{x}\times n_{y}$	Erro ( $\phi$ )	Ordem	Erro ( $u,v$ )	Ordem	Erro ( $\phi$ )	Ordem	Erro ( $u,v$ )	Ordem
$20\times 20$	2.595 $\times 10^{-3}$	–	1.283 $\times 10^{-2}$	–	1.797 $\times 10^{-3}$	–	7.478 $\times 10^{-3}$	–
$40\times 40$	4.231 $\times 10^{-4}$	2.617	2.032 $\times 10^{-3}$	2.658	3.377 $\times 10^{-4}$	2.412	1.169 $\times 10^{-3}$	2.678
$80\times 80$	7.519 $\times 10^{-5}$	2.492	2.924 $\times 10^{-4}$	2.797	6.111 $\times 10^{-5}$	2.466	1.653 $\times 10^{-4}$	2.822
$160\times 160$	1.197 $\times 10^{-5}$	2.651	3.938 $\times 10^{-5}$	2.892	9.602 $\times 10^{-6}$	2.670	2.217 $\times 10^{-5}$	2.898
$320\times 320$	1.735 $\times 10^{-6}$	2.786	5.128 $\times 10^{-6}$	2.941	1.397 $\times 10^{-6}$	2.781	2.895 $\times 10^{-6}$	2.937
$640\times 640$	2.364 $\times 10^{-7}$	2.876	6.554 $\times 10^{-7}$	2.968	1.939 $\times 10^{-7}$	2.849	3.723 $\times 10^{-7}$	2.959

Exemplo 3.2.2: Este problema hipotético bidimensional é um exemplo utilizado por [23,45,46]. Neste problema as velocidades iniciais são todas nulas, a profundidade a montante é de ${\textstyle 10}$ m, enquanto que a profundidade a jusante é assumida ser ${\textstyle 5}$ m ou ${\textstyle 0.1}$ m. O domínio computacional consiste de uma região de ${\textstyle 200{\mbox{m}}\times 200{\mbox{m}}}$ , com uma parede que se estende paralelamente ao eixo ${\textstyle y}$ , tendo ${\textstyle 10}$ m de largura e está centrada em ${\textstyle x=100}$ m. A falha é suposta ser instantânea, possui ${\textstyle 75}$ m de extensão a partir de ${\textstyle y=95}$ m. O canal é horizontal e desconsidera-se resistência ao escoamento. Espera-se a formação de uma frente de choque após o rompimento. Os resultados em ${\textstyle t=7.2}$ s são computados a partir de uma malha uniforme composta por ${\textstyle 101\times 101}$ células e ${\textstyle CFL=0.4}$ . Constam nas Figuras 4a e 4b gráficos da profundidade, vetores velocidade e curvas de nível para o problema com profundidade inicial a jusante de ${\textstyle 5}$ m, enquanto que as Figuras 4c e 4d referem-se ao problema com profundidade inicial a jusante de ${\textstyle 0.1}$ m. Os resultados são consistentes com aqueles presentes na literatura [23,45,46].


(a)	(b)

(c)	(d)
Figura 4. Elevação da superfície da água $h$ , (a) e (c), vetores velocidade e curvas de nível, (b) e (d), para a solução calculada no tempo $t=7.2{\mbox{s}}$ com $CFL=0.4$ , do problema de ruptura de barragem anti-simétrica em um domínio horizontal e sem fricção, com profundidades a montante de $10$ m e a jusante de $5$ m (a)-(b) e $0.1$ m (c)-(d).

Exemplo 3.2.3: Este problema teste visa avaliar a habilidade do esquema em preservar simetria. Considera-se um domínio de ${\textstyle 50{\mbox{m}}\times 50{\mbox{m}}}$ com condições iniciais:

(33)

No instante da falha da barragem, supõe-se que a parede circular seja removida completamente e, subsequentemente, formam-se ondas que se espalham radialmente. A solução numérica é computada numa malha retangular uniforme com ${\textstyle 101\times 101}$ células e o passo de tempo é tal que ${\textstyle CFL=0.6}$ . A Fig. 5a apresenta o perfil da superfície da água ${\textstyle 0.72}$ s após a hipotética falha na barragem circular. Vetores velocidade e curvas nível relativas à superfície ${\textstyle h}$ estão dispostos na Fig. 5b. A simetria da solução numérica é bem preservada e está de acordo com a literatura [7,47].


(a)	(b)
Figura 5. Perfil da superfície da água $h$ (a), vetores velocidade e curvas de nível (b), para a solução calculada no tempo $t=0.72{\mbox{s}}$ com $CFL=0.6$ , do problema de ruptura de barragem circular em um domínio horizontal e sem fricção.

4. Conclusões

Este artigo apresentou o desenvolvimento de um novo esquema CE/SE explícito para a solução das equações de águas rasas em uma e duas dimensões. Na formulação construída, as variáveis dinâmicas e, consequentemente, as leis diferencial e integral de conservação, foram aproximadas localmente por expansões de Taylor de segunda ordem. O conjunto de variáveis de marcha tornou-se constituído pelas variáveis fluxo ${\textstyle q_{m}}$ e suas derivadas espaciais de primeira ${\textstyle (q_{m})_{x}}$ , ${\textstyle (q_{m})_{y}}$ e segunda ordem ${\textstyle (q_{m})_{xx},(q_{m})_{xy}}$ e ${\textstyle (q_{m})_{yy}}$ . Os experimentos computacionais realizados indicaram acurácia de terceira ordem sobre os problemas hiperbólicos lineares testados. As soluções numéricas dos problemas teste de ruptura de barragem unidimensional e bidimensional, caracterizados pela formação de choque, descontinuidade ou simetria, apresentaram concordância com as soluções analíticas e/ou com a literatura. Conclui-se, com isso, que o esquema proposto possui considerável habilidade em capturar choques e descontinuidades fazendo com que seja uma boa ferramenta para análise de fluxo de ruptura de barragem.

Agradecimentos

Os autores agradecem à CAPES pelo apoio financeiro, ao Programa de Pós-Graduação em Métodos Numéricos em Engenharia (PPGMNE - UFPR) e ao Instituto Federal de Educação, Ciência e Tecnologia Catarinense, Campus Araquari, pelo apoio à pesquisa.

Referências

[1] S. Yang, A. Kurganov, Y. Liu, Well-balanced central schemes on overlapping cells with constant subtraction techniques for the Saint-Venant shallow water system. Journal of Scientific Computing 63 (3) (2015) 678–698.

[2] J. Singh, M. S. Altinakar, Y. Ding, Two-dimensional numerical modeling of dam-break flows over natural terrain using a central explicit scheme. Advances in Water Resources 34 (10) (2011) 1366–1375.

[3] J. Soler, E. Bladé, M. Sánchez-Juny, Ensayo comparativo entre modelos unidimensionales y bidimensionales en la modelización de la rotura de una balsa de materiales sueltos erosionables. Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería 28 (2) (2012) 103–111.

[4] H. Ozmen-Cagatay, S. Kocaman, H. Guzel, Investigation of dam-break flood waves in a dry channel with a hump. Journal of Hydroenvironment Research 8 (3) (2014) 304–315.

[5] S. Kocaman, H. Ozmen-Cagatay, Investigation of dam-break induced shock waves impact on a vertical wall, Journal of Hydrology 525 (2015) 1–12.

[6] T.-J. Chang, H.-M. Kao, K.-H. Chang, M.-H. Hsu, Numerical simulation of shallow-water dam break flows in open channels using smoothed particle hydrodynamics. Journal of Hydrology 408 (1-2) (2011) 78–90.

[7] R. Akoh, S. Li, F. Xiao, A CIP/multi-moment finite volume method for shallow water equations with source terms. Int. J. Numer. Meth. Fluids 56 (2008) 2245–2270.

[8] K. Toda, Y. Ogata, T. Yabe, Multi-dimensional conservative semi-Lagrangian method of characteristics CIP for the shallow water equations. Journal of Computational Physics 228 (13) (2009) 4917–4944.

[9] R. Borsche, A well-balanced solver for the Saint Venant equations with variable cross-section. Journal of Numerical Mathematics 23 (2) (2015) 99–115.

[10] J.-W.Wang, R.-X. Liu, Combined finite volume-finite element method for shallow water equations. Computers & Fluids 34 (10) (2005) 1199–1222.

[11] M. C. Díaz, J. López-García, C. Parés, High order exactly well-balanced numerical methods for shallow water systems. Journal of Computational Physics 246 (2013) 242–264.

[12] F. Benkhaldoun, S. Sari, M. Seaid, Projection finite volume method for shallow water flows. Mathematics and Computers in Simulation 118 (2015) 87 – 101, MAMERN V-2013: 5th International Conference on Approximation Methods and Numerical Modeling in Environment and Natural Resources.

[13] S.-J. Liang, T.-W. Hsu, Least-squares finite-element method for shallow-water equations with source terms. Acta Mechanica Sinica 25 (5) (2009) 597–610.

[14] C. Dawson, J. H. Videman, A streamline diffusion finite element method for the viscous shallow water equations. Journal of Computational and Applied Mathematics 251 (2013) 1–7.

[15] G. Kesserwani, D. Caviedes-Voullième, N. Gerhard, S. M Multiwavelet discontinuous Galerkin h-adaptive shallow water model. Computer Methods in Applied Mechanics and Engineering 294 (2015) 56–71.

[16] D. Wirasaet, S. Brus, C. Michoski, E. Kubatko, J. Westerink, C. Dawson, Artificial boundary layers in discontinuous Galerkin solutions to shallow water equations in channels. Journal of Computational Physics 299 (2015) 597–612.

[17] M. Rasulov, Z. Aslan, O. Pakdil, Finite differences method for shallow water equations in a class of discontinuous functions. Applied Mathematics and Computation 160 (2) (2005) 343–353.

[18] X. Zhengfu, S. Chi-Wang, Anti-diffusive finite difference weno methods for shallow water with transport of pollutant. Journal of Computational Mathematics 24 (3) (2006) 239–251.

[19] Y.-H. Hwang, A characteristic particle method for the Saint Venant equations. Computers & Fluids 76 (2013) 58–72.

[20] H. Liu, H. Wang, S. Liu, C. Hu, Y. Ding, J. Zhang, Lattice Boltzmann method for the Saint-Venant equations. Journal of Hydrology 524 (2015) 411–416.

[21] S.-C. Chang, W.-W. To, A new numerical framework for solving conservation laws - the method of space-time conservation element and solution element. Technical Memo TM 104495, NASA, Lewis Research Center, NASA (August 1991).

[22] T. Molls, F. Molls, Space-time conservation method applied to Saint Venant equations. Journal of Hydraulic Engineering 124 (5) (1998) 501–508.

[23] Y. Zhang, Z. Zeng, J. Chen, The improved space-time conservation element and solution element scheme for two-dimensional dam-break flow simulation. International Journal for Numerical Methods in Fluids 68 (5) (2012) 605–624.

[24] S. Qamar, S. Zia, W. Ashraf, The space-time CE/SE method for solving single and two-phase shallow flow models. Computers & Fluids 96 (2014) 136–151.

[25] S.-C. Chang, New developments in the method of space-time conservation element and solution element - aplications to the Euler and Navier-Stokes equations. Technical Memo TM 106226, NASA, Lewis Research Center, NASA (August 1993).

[26] S.-C. Chang, The method of space-time conservation element and solution elemen - a new approach for solving the Navier-Stokes and Euler equations. Journal of Computational Physics 119 (2) (1995) 295–324.

[27] S.-C. Chang, X.-Y. Wang, C.-Y. Chow, The space-time conservation element and solution element method: A new high-resolution and genuinely multidimensional paradigm for solving conservation laws. Journal of Computational Physics 156 (1) (1999) 89–136.

[28] S.-C. Chang, X.-Y. Wang, W.-M. To, Application of the space-time conservation element and solution element method to one-dimensional convection-diffusion problems. Journal of Computational Physics 165 (1) (2000) 189–215.

[29] Z.-C. Zhang, S. T. John Yu, X.-Y.Wang, S.-C. Chang, A. Himansu, P. C. E. Jorgenson, The CE/SE method for Navier-Stokes equations using unstructured meshes for flows at all speeds. AIAA 2000-0393 (2000) 7.

[30] Z.-C. Zhang, S. T. J. Yu, S.-C. Chang, A space-time conservation element and solution element method for solving the two- and threedimensional unsteady Euler equations using quadrilateral and hexahedral meshes. J. Comput. Phys. 175 (1) (2002) 168–199.

[31] S.-C. Chang, A new approach for constructing highly stable high order CESE schemes. Technical Memo TM 2010-216766, NASA, Lewis Research Center, NASA (August 2010).

[32] D. L. Bilyeu, S.-T. J. Yu, Y. Y. Chen, J. L. Cambier, A two-dimensional fourth-order unstrutured-meshed Euler solver based on the CESE method. J. Comput. Phys. 257 (1) (2014) 981–999.

[33] S. Jerez, J. Romero, M. Roselló, F. Arnau, A semi-implicit space-time CE-SE method to improve mass conservation through tapered ducts in internal combustion engines. Mathematical and Computer Modelling 40 (9-10) (2004) 941–951.

[34] S. Qamar, G. Warnecke, Application of space-time CE/SE method to shallow water magnetohydrodynamic equations. Journal of Computational and Applied Mathematics 196 (1) (2006) 132–149.

[35] X. Luo, M. Wang, J. Yang, G. Wang, The space-time CESE method applied to phase transition of water vapor in compressible flows. Computers & Fluids 36 (7) (2007) 1247–1258.

[36] Y. Chou, R.-J. Yang, Application of CESE method to simulate non-Fourier heat conduction in finite medium with pulse surface heating. International Journal of Heat and Mass Transfer 51 (13-14) (2008) 3525–3534.

[37] Z.-C. Zhang, S. T. John Yu, S.-C. Chang, A. Himansu, P. C. E. Jorgenson, A modified space-time conservation element and solution element method for Euler and Navier-Stokes equations. AIAA 99-3277 (1999) 10.

[38] Y.-I. Lim, S.-C. Chang, S. B. Jørgensen, A novel partial differential algebraic equation (PDAE) solver: iterative space-time conservation element/solution element (CE/SE) method. Computers & Chemical Engineering 28 (8) (2004) 1309–1324.

[39] S.-T. J. Yu, L. Yang, R. L. Lowe, S. E. Bechtel, Numerical simulation of linear and nonlinear waves in hypoelastic solids by the CESE method. Wave Motion 47 (3) (2010) 168–182.

[40] S.-C. Chang, X.-Y. Wang, C.-Y. Chow, The space-time conservation element and solution element method - a new high-resolution and genuinely multidimensional paradigm for solving conservaton laws. Technical Memo TM 1998-208843, NASA, Lewis Research Center, NASA (August 1998).

[41] M. Lukácová, K. W. Morton, G. Warnecke, Evolution Galerkin methods for hyperbolic systems in two space dimensions, Mathematics of Computation 69 (232) (2000) 1355–1384.

[42] Qurrat-Ul-Ain, S. Qamar, G. Warnecke, A high-resolution space-time conservative method for non-linear hyperbolic conservation laws. International Journal of Computational Methods 4 (2) (2007) 223–247.

[43] J. J. Stoker, Water Waves: The Mathematical Theory with Applications, Interscience Publishers, New York, 1957.

[44] O. Delestre, C. Lucas, P.-A. Ksinant, F. Darboux, C. Laguerre, T.-N.-T. Vo, F. James, S. Cordier, SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies. International Journal for Numerical Methods in Fluids 72 (3) (2013) 269–300.

[45] R. J. Fennema, M. H. Chaudhry, Implicit methods for two-dimensional unsteady free-surface flows. Journal of Hydraulic Research 27 (3) (1989) 321–332.

[46] C. Zoppou, S. Roberts, Numerical solution of the two-dimensional unsteady dam break. Applied Mathematical Modelling 24 (7) (2000) 457–475.

[47] F. Alcrudo, P. Garcia-Navarro, A high-resolution Godunov-type scheme in finite volumes for 2D shallow-water equations. International Journal for Numerical Methods in Fluids 16 (1993) 489–505.

Latest revision as of 09:58, 10 January 2018

Resumo

Abstract

1 Introdução