The 2D Laplace Equation

Problem Description

$$\frac{\partial^2{p}}{\partial^2{x}}+\frac{\partial^2{p}}{\partial^2{y}} = 0$$

• Analytical Solution:

$$p(x,y) = \frac{x}{4}-4\sum_{m=1(odd)}^{\inf}\frac{\sin{h(m\pi{x})}\cos{(m\pi{y})}}{(m\pi)^2\sin{h(2m\pi{y})}}$$

Problem Formulation

Input Data:

• niter = 51 - maximum number of iterations
• nx = 21 - number of x spatial points
• ny = 11 - number of y spatial points
• xmax = 2
• ymax = 1
• infinity = 100 - a large number

ICs

• $\forall{(x,y)}$ if $n=0 \rightarrow p =0$

BCs

• Dirichlet:

  • $\forall{(n,y)}$ if $x=0 \rightarrow p=0$
  • $\forall{(n,y)}$ if $x=2 \rightarrow p=y$

• Neumann:

  • $\forall{(n,x)}$ if $y=0 \rightarrow \frac{\partial{p}}{\partial{y}}=0$
  • $\forall{(n,x)}$ if $y=1 \rightarrow \frac{\partial{p}}{\partial{y}}=0$

Output Data:

• $\forall{(x, y, n)}p = ?$

Algorthm Design

• $i \rightarrow$ index of lattice in x
• $j \rightarrow$ index of lattice in y
• $n \rightarrow$ index of iterations
• $m \rightarrow$ index of oddd numbers in analytical solution

Algorithm Implementation

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

def laplace_equation_numerical(niter, nx, ny, xmax, ymax):

    """
    Returns the velocity field and distance for 2D linear convection
    """
    # Time steps in space:
    dx = xmax/(nx-1)
    dy = ymax/(ny-1)

    # Data structure initialization:
    p = np.zeros(((nx,ny,niter)))
    x = np.zeros(nx)
    y = np.zeros(ny)
    jpos = np.zeros(ny)
    jneg = np.zeros(ny)

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

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

    # Initial conditions
    p[:,:,0] = 0

    # Dirichlet Boundary Conditions:

    # p boundary, left:
    p[0,:,:] = 0

    # p boundary, right:
    for j in range(0,ny):
        p[nx-1,j,:] = y[j] 

    # von Neumann Boundary Conditions:

    # Values for correction at boundaries:
    for j in range(0,ny-1):
        jpos[j] = j + 1
        jneg[j] = j - 1

    # Set Reflection: 
    jpos[ny-1] = ny-2 
    jneg[0] = 1 

    # Loop
    for n in range(0,niter-1):
        for i in range(1,nx-1):
            for j in range(0,ny):
                p[i,j,n+1] = (( dy**2 * ( p[i+1,j,n]+p[i-1,j,n] ) 
                       + dx**2 * ( p[i,int(jpos[j]),n]+p[i,int(jneg[j]),n] ) ) 
                       / (2*(dx**2 + dy**2)))
            
    return p, x, y

p1, x1, y1 = laplace_equation_numerical(100, 51, 51, 2.0, 1.0)

def plot_2D(p,x,nt,title):
    """
    Plots the 1D velocity field
    """
    plt.figure()
    ax=plt.subplot(111)
    colour=iter(cm.rainbow(np.linspace(0,20,nt)))   
    for n in range(0,nt,20):
        c=next(colour)
        ax.plot(x,p[:,-1,n],linestyle='-',c=c,label='n='+str(n)+' numerical')
    box=ax.get_position()
    ax.set_position([box.x0, box.y0, box.width*1.5,box.height*1.5])
    ax.legend( bbox_to_anchor=(1.02,1), loc=2)
    plt.xlabel('x (m)')
    plt.ylabel('p (Pa)')
    plt.ylim([0,1.0])
    plt.xlim([0,2.0])
    plt.title(title)
    plt.show()
    
def plot_3D(p,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('p (Pa)')
    ax.set_xlim(0,2)
    ax.set_ylim(0,1)
    ax.view_init(30,225)
    Y,X=np.meshgrid(y,x) #note meshgrid uses y,x not x,y!!!
    surf=ax.plot_surface(X,Y,p[:,:,time], rstride=1, cstride=1, cmap=cm.viridis,
            linewidth=0, antialiased=False)
    plt.title(title)
    plt.show()

plot_3D(p1,x1,y1,0,'Figure 1: Initial Conditions: n=0, niter=100, nx=51, ny=51','p (Pa)')
plot_3D(p1,x1,y1,99,'Figure 2: Final Conditions: n=99, niter=100, nx=51, ny=51','p (Pa)')
plot_2D(p1,x1,99,'Figure 3: At Far Neumann Boundary: n=99, niter=100, nx=51, ny=51')