In Module 4, we solved a two-dimentional heat diffusion equation which included Dirichlet and Neumann boundary conditions using an implicit scheme. Now we will solve a similar problem using robin bounday condition.
Robin boundary conditions are a weighted combination of Dirichlet boundary condition and Neumann boundary conditions. This contrasts to the mixed boundary condition, which are boundary conditions of different types specified on different subsets of the boundary. Robin boundary condition are also called impipedance boundary conditions, from their applications in electromagnetic problems, or convective boundary conditions, from their applications in heat transfer problems [1]. $$ k \frac{\partial T}{\partial t}+ hT = hT_{\infty} \quad(1) $$
Lets define our problem. we are using the same grid defined in notebook on 2D Implicit subject to dirichlet boundary condition on the three sides and robin boundary condition on one side. In this note book we will see how the temperature is varied if we cool the silicon chip on one side with air. To refreash, lets rewrite the problem statement with this conditions.
We will consider the chip as a $2D$ plate of size $1{\rm cm}\times 1{\rm cm}$, made of silicon: $\kappa = 159{\rm W/m C}$,$c_p = 0.712\cdot 10^3 {\rm J/kg C}$, $\rho = 2329{\rm kg/m}^3$, and diffusivity $\alpha \approx 10^{-4}{\rm m}^2{/rm s}$. Just to demonstrate a numerical solution. Lets say that the one of the sides of the chip is cooled by air which is has heat transfer coefficent of $h = 32 {\rm W/m^2 k}$, and the velocity of air is $v = 10 {\rm m/s}$.The other three edges are touching other components that have a constant temperature of of $T = 373{\rm K}$ when the machine is operating.The room temperature is $T_{\infty} = 293 {\rm K}$. Now lets see how the temperature varies across the chip as it is cooled from one side.
In [1]:
from IPython.display import display
from IPython.display import Image
display(Image(filename="../heat transfer with robin boundary condition/2dchip (2).jpg"))
Figure 1: Simplified microchip problem setup.
lets begin!!
In Note book on 2D Explicit, we have looked at the heat conduction equation. Lets rewrite the equation again in this note book: $$ \frac{\partial T}{\partial t} = \alpha \left(\frac{\partial^2 T}{\partial x^2}+ \frac{\partial^2 T}{\partial y^2} \right) \quad(2) $$
The implicit scheme for 2D heat equation with central difference in space is written as:
$$ \frac{T^{n+1}_{i,j} - T^n_{i,j}}{\Delta t} = \alpha \left (\frac{T^{n+1}_{i+1, j} - 2T^{n+1}_{i,j} + T^{n+1}_{i-1,j}}{\Delta x^2} + \frac{T^{n+1}_{i, j+1} - 2T^2{n+1}_{i,j} + T^{n+1}_{i,j-1}}{\Delta y^2}\right) \quad(3) $$Rearrange the equation such that all the unknown terms are on the left side and all the known terms are on the right side. And also lets assume that the mesh spacing is equal in both the directions i.e., $\Delta X = \Delta Y = \delta$. $$ -T^{n+1}_{i-1,j} - T^{n+1}_{i+1,j} + \left(\frac{\delta^2}{\alpha \Delta t} + 4 \right) T^{n+1}_{i,j} - T^{n+1}_{i,j-1} - T^{n+1}_{i,j+1} = \frac{\delta^2}{\alpha \Delta t}T^n_{i,j} \quad(4) $$
you can review notebook on 2D Implicit, to understand how to form the grid and how to form the matrix. Since we have learned them in that module, we will directly jump to solving the boundary conditions and substituting them in our 2D heat equation.
In our problem the top, left and bottom boundaries have Dirichlet boundary condition and rignt boundary has robin boundary condition.
Now let's looks at each case,
Bottom boundary condition:
The equation for $j=1$ (interior point adjacent to the bottom boundary) uses values from $j=0$, which are known. Let's put that on the right hand side of the equation. We get this quation for all points across the $x$-axis that are adjacent to teh bottom boundary: $$ -T^{n+1}_{i-1,1} - T^{n+1}_{i+1,1} + \left(\frac{\delta^2}{\alpha \Delta t} + 4 \right) T^{n+1}_{i,1}-T^{n+1}_{i,j+1} = \frac{\delta^2}{\alpha \Delta t}T^n_{i,1} + T^{n+1}_{i,0} \quad(5) $$
Left boundary condition:
similar to bottom boundary condition, the equation for $i=1$ (interior points adjacent to the left boundary) uses known values from $i=0$, and we will put that on the right-hand side:
$$ -T^{n+1}_{2,j} + \left(\frac{\delta^2}{\alpha \Delta t} + 4 \right) T^{n+1}_{1,j} - T^{n+1}_{1,j-1}-T^{n+1}_{1,j+1} = \frac{\delta^2}{\alpha \Delta t}T^n{1,j} + T^{n+1}_{0,j} \quad(6) $$Right boundary condition:
The boundary condition for the right is Robin boundary condition: $$ k \frac{\partial T}{\partial t}+ hT = hT_{\infty} $$ Its finite difference approximation is $$ \frac{T^{n+1}_{n_x-1,j} - T^{n+1}_{n_x-2,j}}{\delta} +hT^{n+1}_{n_x-1,j} = hT_{\infty} \quad(7) $$
Rearranging the terms , we have
$$ T^{n+1}_{n_x-1,j} = \frac{{\delta}hT{\infty}}{k+{\delta}h} +\frac{k}{k+{\delta}h} T^{n+1}_{n_x-2,j} \quad(8) $$for the sake of simplicity lets take $A1 = \frac{{\delta}hT{\infty}}{k+{\delta}h}$ and $B = \frac{k}{k+{\delta}h}$
Finite difference equation for $i = n_x-2$: $$ -T^{n+1}_{n_x-3,j} + \left(\frac{\delta^2}{\alpha \Delta t} + 4 - B \right) T^{n+1}_{n_x-2,j} - T^{n+1}_{n_x-2,j-1}-T^{n+1}_{n_x-2,j+1} = \frac{\delta^2}{\alpha \Delta t}T^n{n_x-2,j} + A1 \quad(9) $$
Top boundary conditions:
Top boundary condition in our case is Dirichlet boundary condition. Therefore the equation for $j = n_y -2 $ we have: $$ -T^{n+1}_{i+1,n_y-2} + \left(\frac{\delta^2}{\alpha \Delta t} + 4 \right) T^{n+1}_{i,n_y-2} - T^{n+1}_{i,n_y-3}-T^{n+1}_{i,n_y-1} = \frac{\delta^2}{\alpha \Delta t}T^n{i,n_y-2} + T^{n+1}_{i-1,n_y-2} \quad(10) $$
Bottom-Left corner:
At $T_{1,1}$ there is Dirichlet boundary condition at $i = 0$ and $j = 0$. The equation for this boundary is as follows:
$$ \small -T^{n+1}_{2,1} + \left(\frac{\delta^2}{\alpha \Delta t} +4 \right) T^{n+1}_{1,1} - T^{n+1}_{1,2} = \frac{\delta^2}{\alpha \Delta t}T^n_{1,1} + T^{n+1}_{0,1} + T^{n+1}_{1,0} \quad(11) $$Top-Left Corner:
At top left corner we have Dirichlet boundary condition at $i = 0$ and $j = n_y-1$. Therefore the equation at point $T_{1,n_y-2}$ is :
$$ \small -T^{n+1}_{2,n_y-2} + \left(\frac{\delta^2}{\alpha \Delta t} +4 \right) T^{n+1}_{2,n_y-2} - T^{n+1}_{1,n_y-3} = \frac{\delta^2}{\alpha \Delta t}T^n_{1,1} + T^{n+1}_{0,ny-2} + T^{n+1}_{1,n_y-1} \quad(12) $$Top-Right Corner:
we have Robin boundary condition at $i = n_x-1$ and Dirichlet boundary condition at $j = n_y-1$. Now at $T_{n_x-2,n_y-2}$ the equation is as follows:
$ \small -T^{n+1}_{n_x-3,n_y-2} + \left(\frac{\delta^2}{\alpha \Delta t} +4 -B \right) T^{n+1}_{n_x-2,n_y-2} - T^{n+1}_{n_x-2,n_y-3}- T^{n+1}_{n_x-2,n_y-1} $ $$ \small \quad = \frac{\delta^2}{\alpha \Delta t}T^n_{n_x-2,n_y-2} + A1 \quad(13) $$
Bottom-Right Corner:
We have Dirichlet boundary condition at the bottom and Robin boundary condition at the top. Therefore the equation for $T_{n_x-2,1}$ is:
$$ \small -T^{n+1}_{n_x-3,1} + \left(\frac{\delta^2}{\alpha \Delta t} +4 -B \right) T^{n+1}_{n_x-2,1} - T^{n+1}_{n_x-2,0}- T^{n+1}_{n_x-2,2} = \frac{\delta^2}{\alpha \Delta t}T^n_{n_x-2,1} +A1 \quad(14) $$we will be solving linear system at every time step.
$$[A][T^{n+1}_\text{int}] = [b]+[b]_{b.c.}$$Since our methodology is same as that in the notebook notebook on 2D Implicit, Figuring out the coefficients for the matrix should be easy.
Now lets start writing the code:
we will be importing scipt.linalg.solve because we need to solve linear system of equations.
In [2]:
import numpy
from scipy.linalg import solve
from math import *
import matplotlib.pyplot as plt
from IPython.display import display
from IPython.display import Image
In [3]:
def constructMatrix(nx, ny, sigma, A1, B):
"""Generate implicit matrix for 2D heat equation with Dirichlet in top, left and bottom and Robin in right
Assume dx = dy
Parameters:
nx : int
number of discrretization points in x
ny : int
number of discretization points in y
sigma: float
alpha*dt/dx
A1 : float
constant in Robin boundary condition
(delta*h*T_out)/((delta*h)+k)
B : float
constant in Robin boundary conditions
(k)/((delta*h)+k)
Return:
-------
A:2D array of floats
Matrix of implicit 2D heat equation
"""
A = numpy.zeros(((nx-2)*(ny-2),(nx-2)*(ny-2)))
row_number = 0 # row counter
for j in range (1,ny-1):
for i in range (1,nx-1):
#corners
if i==1 and j==1:# Bottom left corner (Dirichlet down and left)
A[row_number,row_number] = 1/sigma+4 #set diagnal
A[row_number,row_number+1] = -1 #fetch i+1
A[row_number,row_number+nx-2] = -1 #fetch j+1
elif i==nx-2 and j==1: #Bottom right corner (Dirichlet down,Robin right )
A[row_number,row_number] = 1/sigma+4-B #set diagnal
A[row_number,row_number-1] = -1 #fetch i-1
A[row_number,row_number+nx-2] = -1 #fetch j+1
elif i==1 and j==ny-2: #Top Left corner (Dirichlet up and left)
A[row_number,row_number] = 1/sigma+4 #set diagnal
A[row_number,row_number+1] = -1 #fetch i+1
A[row_number,row_number-(nx-2)] = -1 #fetch j-1
elif i==nx-2 and j==ny-2: #Top right corner (Dirichlet up,Robin right )
A[row_number,row_number] = 1/sigma+4-B #set diagnal
A[row_number,row_number-1] = -1 #fetch i-1
A[row_number,row_number-(nx-2)] = -1 #fetch j-1
# Sides
elif i==1: #Left boundary (Dirichlet)
A[row_number,row_number] = 1/sigma +4 #set diagonal
A[row_number,row_number+1] = -1 #fetch i+1
A[row_number,row_number+nx-2] = -1 #fetch j+1
A[row_number,row_number-(nx-2)] = -1 #fetch j-1
elif i==nx-2: # Right boundary (Robin)
A[row_number,row_number] = 1/sigma +4 - B #set diagonal
A[row_number,row_number-1] = -1 #fetch i-1
A[row_number,row_number+nx-2] = -1 #fetch j+1
A[row_number,row_number-(nx-2)] = -1 #fetch j-1
elif j==1: #Bottom boundary (Dirichlet)
A[row_number,row_number] = 1/sigma +4 #set diagonal
A[row_number,row_number+1] = -1 #fetch i+1
A[row_number,row_number-1] = -1 #fetch i-1
A[row_number,row_number+nx-2] = -1 #fetch j+1
elif j==ny-2: # Top boundary (Dirichlet)
A[row_number,row_number] = 1/sigma +4 #set diagonal
A[row_number,row_number+1] = -1 #fetch i+1
A[row_number,row_number- 1] = -1 #fetch i-1
A[row_number,row_number-(nx-2)] = -1 #fetch j-1
#interior points
else:
A[row_number,row_number] = 1/sigma +4 #set diagonal
A[row_number,row_number+1] = -1 #fetch i+1
A[row_number,row_number-1] = -1 #fetch i-1
A[row_number,row_number+(nx-2)] = -1 #fetch j+1
A[row_number,row_number-(nx-2)] = -1 #fetch j-1
row_number += 1 #jump to next row of the matrix!
return A
In [4]:
def generateRHS(nx, ny, sigma, T, T_bc, A1, B):
""" Generates right-hand side for 2D implicit heat equation with Dirichlet in Top, Left and bottom and Robin in the right
Assume dx=dy,
parameters:
-----------
nx : int
number of discretization points in x
ny : int
number of discretization points in y
sigma: float
alpha*dt/dx
T : array of float
Temperature in current time step
T_bc : float
Temperature in Dirichlet boundary condition
A1 : float
constant in Robin boundary condition
(delta*h*T_out)/((delta*h)+k)
B : float
constant in Robin boundary conditions
(k)/((delta*h)+k)
Return :
-------
RHS : array of float
Right hand side of the 2D implicit heat eqation
"""
RHS = numpy.zeros((nx-2)*(ny-2))
row_number = 0 #row counter
for j in range(1,ny-1):
for i in range(1,nx-1):
#corners
if i==1 and j==1: #Bottom left corner (Dirichlet down and left)
RHS[row_number] = T[j,i]*1/sigma + 2*T_bc
elif i==nx-2 and j==1: # Bottom right corner (Dirichlet down,Robin right)
RHS[row_number] = T[j,i]*1/sigma + T_bc + A1
elif i==1 and j==ny-2: # Top left corner (Dirichlet up and left)
RHS[row_number] = T[j,i]*1/sigma + 2*T_bc
elif i==nx-2 and j==ny-2: # Top right corner (Dirichlet up and Robin right)
RHS[row_number] = T[j,i]*1/sigma + T_bc + A1
#Sides
elif i==1: # Left boundary (Dirichlet)
RHS[row_number] = T[j,i]*1/sigma + T_bc
elif i==nx-2: # Right boundary (Robin)
RHS[row_number] = T[j,i]*1/sigma + A1
elif j==1: # Bottom boundary (Dirichlet)
RHS[row_number] = T[j,i]*1/sigma + T_bc
elif j==ny-2: # Top boundary (Dirichlet)
RHS[row_number] = T[j,i]*1/sigma + T_bc
# Interior point
else:
RHS[row_number] = T[j,i]*1/sigma
row_number += 1 #jump to next row
return RHS
In [5]:
def map_1Dto2D(nx, ny, T_1D, T_bc):
"""Takes temperatures of solution of linear system, stored in 1D,
and puts them in a 2D array with the BCs
Valid for constant Dirichlet bottom, left and bottom and Robin in the right
parameters:
----------
nx : int
number of nodes in x direction
ny : int
number of nides in y direction
T_1D: ARRAY OF FLOATS
SOLUTION OF LINEAR SYSTEM
T_bc: float
Dirichlet BC
Return:
-------
T:2D array of float
Temperature stored in 2D array with BCs
"""
T = numpy.zeros((ny,nx))
row_number = 0
for j in range(1,ny-1):
for i in range(1,nx-1):
T[j,i] = T_1D[row_number]
row_number += 1
# Dirichlet BC
T[0,:] = T_bc
T[:,0] = T_bc
T[-1,:]= T_bc
#Robin BC
T[:,-1] = A1+B*T[:,-2] # Right side
return T
To advanse in time, we will use the following function
In [6]:
def btcs_2D(T, A, nt, sigma, T_bc, nx, ny, dt, A1, B):
"""Advance diffusion equation in time with backward Euler
Parameters:
----------
T: 2D array of float
initial temperature profile
A = 2D array of float
matrix with discretized diffusion equation
nt: int
number of time steps
sigma: float
alpha*dt/dx^2
T_bc :float
Dirichlet BC temperature
nx : int
number of nodes in x direction
ny : int
number of nides in y direction
dt : float
Time step size
A1 : float
constant in Robin boundary condition
(delta*h*T_out)/((delta*h)+k)
B : float
constant in Robin boundary conditions
(k)/((delta*h)+k)
Return:
-------
T: 2D array of floats
temperture profile after nt time steps
"""
j_mid = (numpy.shape(T)[0])/2
i_mid = (numpy.shape(T)[1])/2
counter=0
error=1.
err=[]
while(error>0.001):
Tn = T.copy()
b = generateRHS(nx, ny, sigma, T, T_bc, A1, B)
#use numpy.linalg.solve
T_interior = solve(A,b)
T = map_1Dto2D(nx, ny, T_interior, T_bc)
error=numpy.sqrt(((T - Tn) ** 2).mean())
err.append(error)
counter+=1
return T,error,counter,numpy.asarray(err)
Lets run the code. To do that we have to give the input conditions.
In [7]:
# Inputs to the current problem.
alpha = 1*(10**(-4))
L = 1*(10**(-2))
H = 1*(10**(-2))
nx = 21
ny = 21
nt =300
dx = L/(nx-1)
dy = H/(ny-1)
x = numpy.linspace(0,L,nx)
y = numpy.linspace(0,L,ny)
T_bc = 373
Ti = numpy.ones((ny,nx))*293
delta = 1*(10**(-2))
T_out = 293 #ambient temperature
k = 150
h = 32
B = (k)/((delta*h)+k) # constants in Robin boundary condition
A1 = (delta*h*T_out)/((delta*h)+k) # constants in Robin boundary condition
In [8]:
sigma = 0.25
A = constructMatrix(nx, ny, sigma, A1, B)
In [9]:
T = Ti.copy()
dt = sigma * min(dx, dy)**2/ alpha
T,error,counter,err = btcs_2D(T, A, nt, sigma, T_bc, nx, ny, dt, A1, B)
Lets plot the error for each iteration.
In [10]:
print error
plotx=numpy.linspace(1,counter,counter)
%matplotlib inline
plt.figure()
plt.plot(plotx,err)
plt.ylim([0,10])
plt.xlabel('Iterations')
plt.ylabel('RMS error')
plt.title('Error v/s Iterations')
plt.show()
#print counter
In [11]:
print counter*dt
Lets plot and see the temperature profile!!
In [12]:
import matplotlib.pyplot as plt
%matplotlib inline
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
my, mx = numpy.meshgrid(y,x)
plt.figure(figsize=(7,7))
plt.contourf(my,mx,T,20);
plt.colorbar()
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Temperature Distribution')
Out[12]:
Great! it works. You can compare this result with the implicit scheme for 2D diffusion problem with neumann boundary condtion and see the difference.
You can always check how temperature varies with time by printing the value of T. Though its not needed for our code, but it will help us understand the effect of convection on the chip.
In [13]:
#print T-273
Lets do some thing interesting!
We are going to remove the heat source (i.e.,we are switching of the machine) and see what happens. Its clear that now all the boundary conditions are robin boundary condition (after removing the heat source), now we want to see how the temperature varies . To understand clearly what we are doing it's better to write down the equations again. Note that this experiment is done on the same chip after doing the above condition.
In [14]:
def constructMatrix_2(nx, ny, sigma, A1, B):
"""Generate implicit matrix for 2D heat equation with Robin boundary condition all sides
Assume dx = dy
Parameters:
nx : int
number of discrretization points in x
ny : int
number of discretization points in y
sigma: float
alpha*dt/dx
A1 : float
constant in Robin boundary condition
(delta*h*T_out)/((delta*h)+k)
B : float
constant in Robin boundary conditions
(k)/((delta*h)+k)
Return:
-------
A:2D array of floats
Matrix of implicit 2D heat equation
"""
A = numpy.zeros(((nx-2)*(ny-2),(nx-2)*(ny-2)))
row_number = 0 # row counter
for j in range (1,ny-1):
for i in range (1,nx-1):
#corners
if i==1 and j==1:# Bottom left corner (Robin down and left)
A[row_number,row_number] = 1/sigma+4-2*B #set diagnal
A[row_number,row_number+1] = -1 #fetch i+1
A[row_number,row_number+nx-2] = -1 #fetch j+1
elif i==nx-2 and j==1: #Bottom right corner (Robin down,Robin right )
A[row_number,row_number] = 1/sigma+4-2*B #set diagnal
A[row_number,row_number-1] = -1 #fetch i-1
A[row_number,row_number+nx-2] = -1 #fetch j+1
elif i==1 and j==ny-2: #Top Left corner (Robin up and left)
A[row_number,row_number] = 1/sigma+4-2*B #set diagnal
A[row_number,row_number+1] = -1 #fetch i+1
A[row_number,row_number-(nx-2)] = -1 #fetch j-1
elif i==nx-2 and j==ny-2: #Top right corner (Robin up,Robin right )
A[row_number,row_number] = 1/sigma+4-2*B #set diagnal
A[row_number,row_number-1] = -1 #fetch i-1
A[row_number,row_number-(nx-2)] = -1 #fetch j-1
# Sides
elif i==1: #Left boundary (Robin)
A[row_number,row_number] = 1/sigma +4-B #set diagonal
A[row_number,row_number+1] = -1 #fetch i+1
A[row_number,row_number+nx-2] = -1 #fetch j+1
A[row_number,row_number-(nx-2)] = -1 #fetch j-1
elif i==nx-2: # Right boundary (Robin)
A[row_number,row_number] = 1/sigma +4 - B #set diagonal
A[row_number,row_number-1] = -1 #fetch i-1
A[row_number,row_number+nx-2] = -1 #fetch j+1
A[row_number,row_number-(nx-2)] = -1 #fetch j-1
elif j==1: #Bottom boundary (Robin)
A[row_number,row_number] = 1/sigma +4-B #set diagonal
A[row_number,row_number+1] = -1 #fetch i+1
A[row_number,row_number-1] = -1 #fetch i-1
A[row_number,row_number+nx-2] = -1 #fetch j+1
elif j==ny-2: # Top boundary (Robin)
A[row_number,row_number] = 1/sigma +4-B #set diagonal
A[row_number,row_number+1] = -1 #fetch i+1
A[row_number,row_number- 1] = -1 #fetch i-1
A[row_number,row_number-(nx-2)] = -1 #fetch j-1
#interior points
else:
A[row_number,row_number] = 1/sigma +4 #set diagonal
A[row_number,row_number+1] = -1 #fetch i+1
A[row_number,row_number-1] = -1 #fetch i-1
A[row_number,row_number+(nx-2)] = -1 #fetch j+1
A[row_number,row_number-(nx-2)] = -1 #fetch j-1
row_number += 1 #jump to next row of the matrix!
return A
Again we have to define the function for Right hand side of the matrix as our boundary condition is changed.
In [15]:
def generateRHS_2(nx, ny, sigma, T, T_bc, A1, B):
""" Generates right-hand side for 2D implicit heat equation with Robin boundary condition on all sides
Assume dx=dy,
parameters:
-----------
nx : int
number of discretization points in x
ny : int
number of discretization points in y
sigma: float
alpha*dt/dx
T : array of float
Temperature in current time step
T_bc : float
Temperature in Dirichlet boundary condition
A1 : float
constant in Robin boundary condition
(delta*h*T_out)/((delta*h)+k)
B : float
constant in Robin boundary conditions
(k)/((delta*h)+k)
Return :
-------
RHS : array of float
Right hand side of the 2D implicit heat eqation
"""
RHS = numpy.zeros((nx-2)*(ny-2))
row_number = 0 #row counter
for j in range(1,ny-1):
for i in range(1,nx-1):
#corners
if i==1 and j==1: #Bottom left corner (Robin down and left)
RHS[row_number] = T[j,i]*1/sigma + 2*A1
elif i==nx-2 and j==1: # Bottom right corner (Robin down,Robin right)
RHS[row_number] = T[j,i]*1/sigma + 2*A1
elif i==1 and j==ny-2: # Top left corner (Robin up and left)
RHS[row_number] = T[j,i]*1/sigma + 2*A1
elif i==nx-2 and j==ny-2: # Top right corner (Robin up and Robin right)
RHS[row_number] = T[j,i]*1/sigma +2* A1
#Sides
elif i==1: # Left boundary (Robin)
RHS[row_number] = T[j,i]*1/sigma + A1
elif i==nx-2: # Right boundary (Robin)
RHS[row_number] = T[j,i]*1/sigma + A1
elif j==1: # Bottom boundary (Robin)
RHS[row_number] = T[j,i]*1/sigma + A1
elif j==ny-2: # Top boundary (Robin)
RHS[row_number] = T[j,i]*1/sigma + A1
# Interior point
else:
RHS[row_number] = T[j,i]*1/sigma
row_number += 1 #jump to next row
return RHS
In [16]:
def map_1Dto2D_2(nx, ny, T_1D, T_bc):
"""Takes temperatures of solution of linear system, stored in 1D,
and puts them in a 2D array with the BCs
Valid for Robin boundary condition on all sides
parameters:
----------
nx : int
number of nodes in x direction
ny : int
number of nides in y direction
T_1D: ARRAY OF FLOATS
SOLUTION OF LINEAR SYSTEM
T_bc: float
Dirichlet BC
Return:
-------
T:2D array of float
Temperature stored in 2D array with BCs
"""
T = numpy.zeros((ny,nx))
row_number = 0
for j in range(1,ny-1):
for i in range(1,nx-1):
T[j,i] = T_1D[row_number]
row_number += 1
#Robin BC
T[0,:] = A1+B*T[1,:]
T[:,0] = A1+B*T[:,1]
T[-1,:]= A1+B*T[-2,:]
T[:,-1] = A1+B*T[:,-2]# Right
return T
In [17]:
def btcs_2D_2(T, A, nt, sigma, T_bc, nx, ny, dt, A1, B):
"""Advance diffusion equation in tim with backward Euler
Parameters:
----------
T: 2D array of float
initial temperature profile
A = 2D array of float
matrix with discretized diffusion equation
nt: int
number of time steps
sigma: float
alpha*dt/dx^2
T_bc :float
Dirichlet BC temperature
nx : int
number of nodes in x direction
ny : int
number of nides in y direction
dt : float
Time step size
A1 : float
constant in Robin boundary condition
(delta*h*T_out)/((delta*h)+k)
B : float
constant in Robin boundary conditions
(k)/((delta*h)+k)
Return:
-------
T: 2D array of floats
temperture profile after nt time steps
"""
j_mid = (numpy.shape(T)[0])/2
i_mid = (numpy.shape(T)[1])/2
Ti = numpy.ones((ny,nx))*293
counter=0
error=1.
err=[]
while(error>0.001):
Tn = T.copy()
b = generateRHS_2(nx, ny, sigma, T, T_bc, A1, B)
#use numpy.linalg.solve
T_interior = solve(A,b)
T = map_1Dto2D_2(nx, ny, T_interior, T_bc)
error=numpy.sqrt(((T - Ti) ** 2).mean())
err.append(error)
counter+=1
return T,error,counter,numpy.asarray(err)
In [18]:
nt =600
A = constructMatrix_2(nx, ny, sigma, A1, B)
dt = sigma * min(dx, dy)**2/ alpha
T,error,counter,err = btcs_2D_2(T, A, nt, sigma, T_bc, nx, ny, dt, A1, B)
In [19]:
print error
plotx=numpy.linspace(1,counter,counter)
%matplotlib inline
plt.figure()
plt.plot(plotx,err)
plt.ylim([0,100])
plt.xlabel('Iterations')
plt.ylabel('RMS error')
plt.title('Error v/s Iterations')
plt.show()
#print counter
In [20]:
#print counter*dt # Time taken for temperature to reach 20 degrees celcius
In [21]:
my, mx = numpy.meshgrid(y,x)
plt.figure(figsize=(7,7))
plt.contourf(my,mx,T,20);
plt.colorbar()
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Temperature Distribution')
Out[21]:
In [22]:
#print T-273
my, mx = numpy.meshgrid(y,x)
plt.figure(figsize=(7,7))
plt.contourf(my,mx,T-Ti,20)
plt.colorbar()
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Error')
Out[22]:
Isn't it cool. If you want to get clear idea, you can always print the value of T and see what happened to the Temperature in our experiment. You must have noticed by now that we have incuded a condition in our code that makes this code run as long as error is greater than 0.001. Which means that, in our experiment as soon as the heat source is removed we are allowing the chip to cool till the temperature reaches close to the atmospheric temperature under free convection. The result show us how the Temperature varies across the chip.
Gustafson, K., (1998). Domain Decomposition, Operator Trignometry, Robin Condition, Contemporary Mathematics, 218. 432-437
Barba, Lorena A., et al. "MEA 6286 Practical Numerical Methods with Python," GW Open edX, The George Washington University, 2014. http://openedx.seas.gwu.edu/courses/GW/MAE6286/2014_fall/about
In [24]:
from IPython.core.display import HTML
css_file = '../../styles/numericalmoocstyle.css'
HTML(open(css_file, "r").read())
Out[24]: