2D Non-Linear Convection Diffusion

$$\frac{\partial{u}}{\partial{t}}+u\frac{\partial{u}}{\partial{x}}++v\frac{\partial{u}}{\partial{y}}=\nu (\frac{\partial^2{u}}{\partial{x^2}}+\frac{\partial^2{u}}{\partial{y^2}})$$$$\frac{\partial{v}}{\partial{t}}+u\frac{\partial{v}}{\partial{x}}+u\frac{\partial{v}}{\partial{y}}=\nu (\frac{\partial^2{v}}{\partial{x^2}}+\frac{\partial^2{v}}{\partial{y^2}})$$

Formulate the Problem

A. Input Data

• $nt$ = 51
• $tmax$ = 0.5
• $nx$ = 21
• $nj$ = 21
• $xmax$ = 2
• $ymax$ = 2
• $nu$ = 0.1

Initial Conditions(ICs) - $t$ = 0

Boundary Conditions (BCs)

• $x$ = 0 and $x$ = 2

• $y$ = 0 and $y$ = 2

• $0\leq{t}\leq0.5 \rightarrow 0\leq{n}\leq50 \rightarrow u(x,y,t) | v(x,y,t) = 1$

B. Output Data

• $0\leq{x}\leq2 \rightarrow 0\leq{y}\leq2 \rightarrow 0\leq{t}\leq0.5$ | $u(x,y,t) | v(x,y,t) = ?$

Implement Algorithm in Python

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.cm as cm

def twod_burgers(nt, nx, ny, tmax, xmax, ymax, nu):
    """
    Returns the velocity field and distance for 2D non-linear convection-diffusion
    """
    # Increments
    dt = tmax/(nt-1)
    dx = xmax/(nx-1)
    dy = ymax/(ny-1)

    # Initialise data structures
    u = np.zeros(((nx,ny,nt)))
    v = np.zeros(((nx,ny,nt)))
    x = np.zeros(nx)
    y = np.zeros(ny)

    # Boundary conditions
    u[0,:,:] = u[nx-1,:,:] = u[:,0,:] = u[:,ny-1,:] = 1
    v[0,:,:] = v[nx-1,:,:] = v[:,0,:] = v[:,ny-1,:] = 1

    # Initial conditions
    u[:,:,:] = v[:,:,:] = 1
    u[int((nx-1)/4):int((nx-1)/2),int((ny-1)/4):int((ny-1)/2),0] = 2
    v[int((nx-1)/4):int((nx-1)/2),int((ny-1)/4):int((ny-1)/2),0] = 2

    # Loop
    for n in range(0,nt-1):
        for i in range(1,nx-1):
            for j in range(1,ny-1):
                u[i,j,n+1] = (u[i,j,n]+
                              -dt*((u[i,j,n]/dx)*(u[i,j,n]-u[i-1,j,n])+
                                   (v[i,j,n]/dy)*(u[i,j,n]-u[i,j-1,n]))+
                              dt*nu*( ( u[i-1,j,n]-2*u[i,j,n]+u[i+1,j,n])/dx**2+
                                     (u[i,j-1,n]-2*u[i,j,n]+u[i,j+1,n])/dy**2))
                v[i,j,n+1] = (v[i,j,n]+
                              -dt*((u[i,j,n]/dx)*(v[i,j,n]-v[i-1,j,n])+
                                   (v[i,j,n]/dy)*(v[i,j,n]-v[i,j-1,n]))+
                              dt*nu*((v[i-1,j,n]-2*v[i,j,n]+v[i+1,j,n])/dx**2+
                                     (v[i,j-1,n]-2*v[i,j,n]+v[i,j+1,n])/dy**2))

    # X Loop
    for i in range(0,nx):
        x[i] = i*dx

    # Y Loop
    for j in range(0,ny):
        y[j] = j*dy

    return u, v, x, y

def plot(u,x,y,time,title,label):
    fig=plt.figure(figsize=(11,7),dpi=100)
    ax=fig.gca(projection='3d')
    ax.set_xlabel('x (m)')
    ax.set_ylabel('y (m)')
    ax.set_zlabel(label)
    X,Y=np.meshgrid(x,y)
    surf=ax.plot_surface(X,Y,u[:,:,time], cmap=cm.viridis, rstride=2, cstride=2)
    plt.title(title)
    #plt.show()
    plt.savefig('burgers.pdf')

u,v,x,y = twod_burgers(311, 51, 51, 0.5, 2.0, 2.0, 0.1)

plot(u,x,y,0,'Figure 1: nu=0.1, nt=51, nx=21, ny=21, t=0sec','u (m/s)')

plot(u,x,y,310,'Figure 2: nu=0.1, nt=51, nx=21, ny=21, t=0.5sec','u (m/s)')

plot(v,x,y,0,'Figure 3: nu=0.1, nt=51, nx=21, ny=21, t=0sec','v (m/s)')

plot(v,x,y,310,'Figure 4: nu=0.1, nt=51, nx=21, ny=21, t=0.5sec','v (m/s)')