Comparation of numerical and analytical solution of 1-D Heat Equation

Introduction

1-d Linear heat Equation:
u t = c 2 u x 2 \frac{\partial u}{\partial t} = c \frac{\partial^2 u}{\partial x^2}

Initial Condition:
f ( x ) = { 0 0<x<0.5 1 0.5<x<1 0 1<x<2 f(x)= \begin{cases} 0& \text{0<x<0.5}\\ 1& \text{0.5<x<1}\\ 0& \text{1<x<2}\\ \end{cases}

Dirichlet Boundary Condition:
u ( 0 , t ) = 0 u ( 2 , t ) = 0 u(0,t) = 0\\ u(2,t) = 0

Design Algorithm to Solve Problem

Discrete equation: Forward Difference in time, Central Difference in space.
u i n + 1 u i n Δ t = c u i + 1 n 2 u i n + u i 1 n Δ x \frac{u_i^{n+1}-u_{i}^{n}} {\Delta t} = c \frac{u_{i+1}^{n}-2u_i^n+u_{i-1}^{n}} {\Delta x}
u i n + 1 = u i n + c Δ t Δ x 2 ( u i + 1 n 2 u i n + u i 1 n ) u_i^{n+1} = u_i^{n} + c \frac{\Delta t}{\Delta x^2} (u_{i+1}^{n} - 2u_i^{n}+u_{i-1}^{n})

code



def diffusion(nt, nx, tmax, length, c):

   dt = tmax/(nt-1)
   dx = length/(nx-1) 
   import numpy as np 
  u = np.zeros((nx,nt))
  x = np.zeros(nx)
 
   # BC

   u[0,:] = u[nx-1,:]= 0

   # IC

   for i in range(1,nx-1):

      if(i*dx > 0.5 and i*dx < 1):

         u[i,0] = 1

      else:

         u[i,0] = 0


   # Loop

   for n in
range(0,nt-1):

      for i in
range(1,nx-1):

        u[i,n+1] =
u[i,n] + c*(dt/dx**2.0)*(u[i+1,n]-2.0*u[i,n]+u[i-1,n]) 

   # X Loop

   for i in range(0,nx):
      x[i] = i*dx
   return u, x


import matplotlib.pyplot as plt
def plot_diffusion(u,x,nt,title):
plt.figure()
for i in range(0,nt,10):
    plt.plot(x,u[:,i],)
    plt.xlabel('x(m)')
    plt.ylabel('u(m/s)')
    plt.ylim([0,2.2])
    plt.title(title)

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
You may find that the solution does not always converge.
I will explain that later (Hopefully)

Analytical solution:

在这里插入图片描述在这里插入图片描述

def ana_diffusion(nt, nx, tmax, length, c):
   dt = tmax/(nt-1)
   dx = length/(nx-1)
   import numpy as np
   u = np.zeros((nx,nt))
   x = np.zeros(nx)
   # X Loop
   for i in range(0,nx):
      x[i] = i*dx
# Loop
   for n in range(0,nt-1):
      for i in range(1,nx-1):
            for k in range(1,1000):
                  u[i,n] = u[i,n] + 2/(k*math.pi)*(math.cos(k*math.pi*0.25)-math.cos(k*math.pi*0.5))*math.exp(-0.1*0.25*k**2*math.pi**2*dt*n)*math.sin(0.5*k*math.pi*x[i])
    return u, x
    

在这里插入图片描述
For analytical solution with Newmann & Robin BC I will do them later (hopefully)

猜你喜欢

转载自blog.csdn.net/Cappucccccino/article/details/106185791