2D Linear Convection

$$\frac{\partial{u}}{\partial{t}}+c\frac{\partial{u}}{\partial{x}}+c\frac{\partial{u}}{\partial{y}}=0$$$$\frac{\partial{v}}{\partial{t}}+c\frac{\partial{v}}{\partial{x}}+c\frac{\partial{v}}{\partial{y}}=0$$

Initial Conditions (ICs):

$$u(x,y) = \begin{cases} \begin{matrix} 2\ \text{for} & 0.5 \leq x, y \leq 1 \cr 1\ \text{for} & \text{everywhere else}\end{matrix}\end{cases}$$

and Boundary Conditions(BCs):

$$u = 1\ \text{for } \begin{cases} \begin{matrix} x = 0,\ 2 \cr y = 0,\ 2 \end{matrix}\end{cases}$$

Implement Algorithm in Python

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

def twodconvection(nt, nx, ny, tmax, xmax, ymax, c):
    """
    Returns the velocity field for 2D convection
    """
    # Increments in space and time
    dt = tmax/(nt-1)
    dx = xmax/(nx-1)
    dy = ymax/(ny-1)
    
    # Initialization of Data Structures
    u = np.ones((nx, ny, nt))
    v = np.ones((nx, ny, nt))
    x = np.zeros(nx)
    y = np.zeros(ny)
    
    # BCs
    u[0,:,:]=u[nx-1,:,:]=u[:,0,:]=u[:,ny-1,:]=1
    v[0,:,:]=v[nx-1,:,:]=v[:,0,:]=v[:,ny-1,:]=1
    
    # ICs
    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]-c*dt*((1/dx)*(u[i,j,n]-u[i-1,j,n])+
                                         (1/dy)*(u[i,j,n]-u[i,j-1,n])))
                v[i,j,n+1] = (v[i,j,n]-c*dt*((1/dx)*(v[i,j,n]-v[i-1,j,n])+
                                         (1/dy)*(v[i,j,n]-v[i,j-1,n])))

    # 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):
    """
    Plots the 2D velocity field
    """
    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],rstride=2,cstride=2)
    plt.title(title)
    #plt.show()
    plt.savefig('filename.pdf')

u,v,x,y = twodconvection(101, 81, 81, 0.5, 2.0, 2.0, 0.5)
plot(u,x,y,0,'Figure 1: c=0.5m/s, nt=101, nx=81, ny=81, t=0sec','u (m/s)')

plot(u,x,y,100,'Figure 2: c=0.5m/s, nt=101, nx=81, ny=81, t=0.5sec','u (m/s)')

Conclusions

• Why isn’t the square wave maintained?

  • As with 1D, the first order backward differencing scheme in space creates false diffusion.

• Why does the wave shift?

  • The square wave is being convected by the constant linear wave speed c in 2 dimensions.
  • For c>0 profiles shift by $c\Delta{t}$ - compare Figure 1, 2, 3 and 4

To save high resolution images:

• save pdf
• and then converting this pdf to a png on the command line so you can useit in powerpoint:

  pdftoppm -png -r 300 filename.pdf filename