Código 4 - Forma Não Conservativa das Equações de Águas Rasas

Definindo os parâmetros iniciais do programa.

# Define a região espacial
Lx = 100.0   
Ly = 100.0   

#define o número de pontos na direção x e y
px = 201     
py = 201     

#define o espaçamento entre os pontos mais próximos
dx = Lx / (px - 1)  
dy = Ly / (py - 1)  

# Localiza o número de pontos ao longo das direções x e y
x = np.linspace(0.0, Lx, num=px)
y = np.linspace(0.0, Ly, num=py)

# Plano de pontos do oceano definidos
X, Y = np.meshgrid(x,y) 

#profundidade do oceano h = 50 m
h = 50 * np.ones_like(X)

# Condição inicial da onda gerado pelo sísmo 
eta0 = 2. * np.exp(-((X-50)**2/10)-((Y-50)**2/10))

# Condição inicial dos fluxos de descarga
M0 = 10. * eta0
N0 = 10. * eta0

# constante da gravidade e de fricção
g = 9.81  
n = 0.025

Tmax = 5.
dt = 0.01
nt = (int)(Tmax/dt)

eta, M, N = shallow_water(eta0,M0,N0,h,g,n,nt,dx,dy,dt,X,Y)

A função shallow water waves recebe os parâmetros iniciais do nosso programa, executa o Loop responsável pela atualização das variáveis da amplitude da onda e do fluxo de descarga com o tempo, através da chamada das funções atualiza M,N e eta. Posteriormente, a cada 10 passagens dentro do loop um plot do sistema é feito.

def shallow_water(eta0, M0, N0, h, g, n, nt, dx, dy, dt, X, Y):
    
    
    # Copia o campos de amplitude da onda, e os fluxos de descarga  
    eta = eta0.copy()
    M = M0.copy()
    N = N0.copy()
    
    D = eta + h
    
    # Retorna o número de pontos na direção x e y para eta
    ny, nx = eta.shape
    
    # Dividi o espaço Lx e Ly definindo os pontos nas direções x e y
    x = np.linspace(0, nx*dx, num=nx)
    y = np.linspace(0, ny*dy, num=ny)
    
    limites = [np.min(x), np.max(x), np.min(y), np.max(y)]
    
    # PLOT
    plt.clf()
    
    
    fig = plt.figure(figsize=(10.,6.))
    plt.tight_layout()
    
    plt.suptitle('Evolução da onda - forma disspitava - n = 0.025', fontsize=14)
    plt.title(f'Tempo: 0 s', fontsize=10)
    
    # plot profundidade
    fundo = plt.imshow(-h, cmap ='Purples', interpolation = 'nearest', extent = limites)
    amp = plt.imshow(eta, extent = limites, interpolation = 'sinc', cmap = 'seismic', alpha= .75, vmin=-1.0, vmax= 1.0)
    
    #plt.plot(f'Tempo 0 s')
    plt.xlabel('x [m]')
    plt.ylabel('y [m]')
    cbar_amp = plt.colorbar(amp)
    cbar_fundo = plt.colorbar(fundo)
    cbar_fundo.set_label(r'$-h$ [m]')
    cbar_amp.set_label(r'$\eta$ [m]')
    #cbar_amp_eta.setlabel = plt.colorbar(r'$eta$ [m]')
    
    plt.show()
    plt.close()
    
    conta = 0
    
    # Loop over timesteps
    for k in range(1,nt):
        

        eta = atualiza_eta(eta, M, N, dx, dy, dt, nx, ny)
        M = atualiza_M(eta, M, N, D, g, h, n, dx, dy, dt, nx, ny)
        N = atualiza_N(eta, M, N, D, g, h, n, dx, dy, dt, nx, ny)

        D = eta + h
        
        conta = conta + 1
        
        if(conta == 10):
            
            plt.clf()
            fig = plt.figure(figsize=(10.,6.))

             # # titulos
            plt.suptitle('Evolução da onda - forma disspitava - n = 0.025', fontsize=14)
            plt.title(f'Tempo: {round(dt*k, 3)} s', fontsize=10)

            fundo = plt.imshow(-h, 'Purples', interpolation = 'nearest', extent = limites)
            amp = plt.imshow(eta, extent = limites, interpolation = 'sinc', cmap = 'seismic', alpha= 0.75, vmin=-1.0, vmax= 1.0)
    
            #plt.title('tempo = %f', dt*n )
            #plt.plot(f'Tempo {round(k*dt,3)} s')
            plt.xlabel('x [m]')
            plt.ylabel('y [m]')
            cbar_amp = plt.colorbar(amp)
            cbar_fundo = plt.colorbar(fundo)
            cbar_fundo.set_label(r'$-h$ [m]')
            cbar_amp.set_label(r'$\eta$ [m]')
    
            conta = 0
            plt.show()
            plt.close()
                
    return eta, M, N

Calculo das equações de águas rasas não conservativa por FTCS.

def atualiza_eta(eta, M, N, dx, dy, dt, nx, ny):
    
    for j in range(1,nx-1):
        for i in range(1,ny-1):
            
            dMdx = (M[j+1,i] - M[j-1,i]) / (2. * dx)
            dNdy = (N[j,i+1] - N[j,i-1]) / (2. * dy)
            eta[j, i] = eta[j, i] - dt * (dMdx + dNdy)
    
    #Condições de contorno do problema
    eta[0,:] = eta[1,:]
    eta[-1,:] = eta[-2,:]
    eta[:,0] = eta[:,1]
    eta[:,-1] = eta[:,-2]
    
    return eta          

def atualiza_M(eta, M, N, D, g, h, n, dx, dy, dt, nx, ny):
    
    M2 = M **2 / D
    MN = M * N / D
    fric = g * n**2 * M * np.sqrt(M**2 + N**2) / D**(7./3.)
    
    for j in range(1,nx-1):
        for i in range(1,ny-1):            
            
            dM2dx = (M2[j+1,i] - M2[j-1,i]) / (2. * dx)
            dMNdy = (MN[j,i+1] - MN[j,i-1]) / (2. * dy)
            dETAdx = (eta[j+1,i] - eta[j-1,i]) / (2. * dx)
            
            M[j, i] = M[j, i] - dt * (dM2dx + dMNdy + g * D[j,i] * dETAdx + fric[j,i])
            
    return M   

def atualiza_N(eta, M, N, D, g, h, n, dx, dy, dt, nx, ny):
    
    
    MN = M * N / D
    N2 = N**2 / D
    fric = g * n**2 * N * np.sqrt(M**2 + N**2) / D**(7./3.)
    
    for j in range(1,nx-1):
        for i in range(1,ny-1):            
            
            dMNdx = (MN[j+1,i] - MN[j-1,i]) / (2. * dx)
            dN2dy = (N2[j,i+1] - N2[j,i-1]) / (2. * dy)
            dETAdy = (eta[j,i+1] - eta[j,i-1]) / (2. * dy)
            
            N[j, i] = N[j, i] - dt * (dMNdx + dN2dy + g * D[j,i] * dETAdy + fric[j,i])
                        
    return N