'1-DECLARAÇÃO DE VARIÁVEIS E ATRIBUIÇÕES INICIAIS: '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Dim Q(300), y(300), z(300), A(300) As Variant
Dim Q1(300), A1(300) As Variant Dim Q2(300), A2(300) As Variant Dim Qo(300), Ao(300) As Variant
Dim H(300), H1(300), b(300), Pm(300) As Variant Dim visc(300) As Variant
So = Sheets("1").Cells(29, 2) L = Sheets("1").Cells(28, 2) nx = 200 dx = L / (nx - 1) zo = 865 base = Sheets("1").Cells(27, 2) k = 3 dt1 = Sheets("1").Cells(19, 2)
td = Sheets("1").Cells(18, 2) * 3600 'duração da chuva em segundos CFL = 0.25
'dt = 1
Sheets("3").Activate
Range(Cells(22, 1), Cells(10000, 7)).ClearContents '2-INICIALIZAÇÃO DAS VARIÁVEIS
'>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> '2.1-Largura inicial do canal: '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> For j = 1 To nx
b(j) = base Next
'2.2-Variáveis dependentes iniciais: equação de Manning '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> For j = 1 To nx z(j) = zo - So * (j - 1) * dx y(j) = 0.001 Q(j) = 0# A(j) = b(j) * y(j) Pm(j) = 2 * y(j) + b(j) H(j) = z(j) + y(j) + (Q(j) / A(j)) ^ 2 / 19.62 Next '3-LAÇO TEMPORAL tempo = 0 nivel = 0 ind = 1 fatemp = 0 '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Do
'3.0-Variáveis no instante atual '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> For j = 1 To nx
Qo(j) = Q(j) Ao(j) = A(j)
y(j) = Ao(j) / b(j)
Next
'Hidrograma de entrada: If (tempo <= ind * dt1) Then If (ind = 1) Then
Tba = 0 Tci = dt1 Qba = 0
Qci = Sheets("1").Cells(ind + 31, 5)
Qent = Qba + (Qci - Qba) * (tempo - Tba) / (Tci - Tba) Else
Tba = (ind - 1) * dt1 Tci = ind * dt1
Qba = Sheets("1").Cells(ind - 1 + 31, 5) Qci = Sheets("1").Cells(ind + 31, 5)
Qent = Qba + (Qci - Qba) * (tempo - Tba) / (Tci - Tba) End If
Qo(1) = Qent End If
If (tempo > ind * dt1) Then
Sheets("3").Cells(ind + 21, 6) = tempo Sheets("3").Cells(ind + 21, 7) = Qo(nx) ind = ind + 1
End If
'3.1-Condição CFL: estabilidade numérica '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> dtmin = 1000 For j = 1 To nx dt = CFL * dx / (Q(j) / b(j) / y(j) + Sqr(9.81 * y(j))) If (dt < dtmin) Then dtmin = dt End If Next dt = dtmin 'dt = 0.01 tempo = tempo + dt
'3.2-ESQUEMA EXPLÍCITO DE MC CORMACK: etapa preditora progressiva '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> For j = 1 To (nx - 1)
'Saint-Venant:
A1(j) = Ao(j) - dt / dx * (Qo(j + 1) - Qo(j))
Q1(j) = Qo(j) - dt / dx * (Qo(j + 1) ^ 2 / Ao(j + 1) - Qo(j) ^ 2 / Ao(j)) - 9.81 * Ao(j) * dt / dx * (y(j + 1) - y(j)) + 9.81 * Ao(j) * dt / dx * (z(j + 1) - z(j)) - 9.81 * Ao(j) * dt / dx * (H(j + 1) - H(j))
H1(j) = z(j) + A1(j) / b(j) + (Q1(j) / A1(j)) ^ 2 / 19.62 y(j) = A1(j) / b(j)
Next
'Condição de contorno na saída: derivadas nulas A1(nx) = A1(nx - 1)
Q1(nx) = Q1(nx - 1) y(nx) = A1(nx) / b(nx)
H1(nx) = z(nx) + A1(nx) / b(j) + (Q1(nx) / A1(nx)) ^ 2 / 19.62 '3.3-ESQUEMA EXPLÍCITO DE MC CORMACK: etapa corretora regressiva '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> For j = 2 To (nx)
'Saint-Venant:
A2(j) = A1(j) - dt / dx * (Q1(j) - Q1(j - 1))
Q2(j) = Q1(j) - dt / dx * (Q1(j) ^ 2 / A1(j) - Q1(j - 1) ^ 2 / A1(j - 1)) - 9.81 * A1(j) * dt / dx * (y(j) - y(j - 1)) + 9.81 * A1(j) * dt / dx * (z(j) - z(j - 1)) - 9.81 * A1(j) * dt / dx * (H1(j) - H1(j - 1))
Next
'Condições de contorno na entrada: A2(1) = 2 * A2(2) - A2(3)
Q2(1) = Qent
'3.4-ESQUEMA EXPLÍCITO DE MC CORMACK: etapa média '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> For j = 1 To (nx)
A(j) = 0.5 * (Ao(j) + A2(j)) Q(j) = 0.5 * (Qo(j) + Q2(j)) y(j) = A(j) / b(j)
H(j) = z(j) + y(j) + (Q(j) / A(j)) ^ 2 / 19.62 Next
'3.4a - CORREÇÃO DE OSCILAÇÕES NUMÉRICAS PELA VISCOSIDADE ARTIFICIAL: '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> For j = 2 To (nx - 1)
visc(j) = Abs(y(j + 1) - 2 * y(j) + y(j - 1)) / Abs(y(j + 1) + 2 * Abs(y(j)) + Abs(y(j - 1)))
maxi = visc(j)
If (visc(j + 1) > maxi) Then maxi = visc(j + 1) End If
e_mais = k * maxi maxi = visc(j)
If (visc(j - 1) > maxi) Then maxi = visc(j - 1) End If
e_menos = k * maxi 'correção:
y(j) = y(j) + e_mais * (y(j + 1) - y(j)) - e_menos * (y(j) - y(j - 1)) Next
y(1) = 2 * y(2) - y(3) y(nx) = y(nx - 1) '3.5-SAÍDA DE DADOS: '>>>>>>>>>>>>>>>>>>>>>>>>> If (tempo > fatemp) Then For i = 1 To nx Cells(21 + i, 1) = (i - 1) * dx Cells(21 + i, 2) = y(i) Cells(21 + i, 3) = z(i) Cells(21 + i, 4) = Q(i) Next fatemp = fatemp + 600 End If Cells(5, 2) = tempo Cells(6, 2) = Q(nx) Cells(7, 2) = Qent tempo = tempo + dt nivel = nivel + 1
Loop Until (Abs(Q(nx)) < 0.2 And tempo > td)
'>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Sheets("1").Activate