I'm trying to write a Python code that is a numerical solver for 1-d heat conduction (using FVM) with a temperature dependent thermal conductivity.
The solver has three functions I need to iterate until convergence. This is basically the algorithm from Patankar's Numerical Heat Transfer and Fluid Flow book.
Conductivity function: takes an initial guess for temperature field and outputs nodal conductivities, dependent on past value of nodal temperatures.
Nodal coefficients: calculates nodal coefficients based on current nodal conductivities.
Nodal temperatures: uses a TDMA algorithm to calculate the nodal temperatures based on current coefficients.
Once the temperatures are calculated, the values are to be cycled through to the conductivity function to continue the loop until convergence.
I can make this work just fine when the conductivity is constant, but start getting issues when I try to implement this pseudocode. I get a maximum recursion depth exceeded error. My code is below:
import numpy as np
import matplotlib.pyplot as plt
import copy
# SOLUTION ALGORITHM
# Step 1: DATA
def data(n, hw, he, TinfW, TinfE):
n = 5
hw = 10**10
he = 10**-10
TinfW = 100
TinfE = 10**12
return n, hw, he, TinfW, TinfE
# Step 2: GRID
def grid(x1, x2, data):
x1 = 0
x2 = 0.5
xe = (x2-x1) / data(n, hw, he, TinfW, TinfE)[0]
return xe
# Step 3: GAMSOR
def gamsor(k, q, TDMA):
k0 = 100
beta = 20
T0 = 25
for i in range(data(n, hw, he, TinfW, TinfE)[0]):
k[i] = k0 + beta#*(TDMA(data, coeff)[i]-T0)
q[i] = 10**5
return k, q
#Step 4: COEFF
def coeff(Ae, Aw, Ap, data, gamsor, grid):
for i in range(data(n, hw, he, TinfW, TinfE)[0]):
Ae[i] = Aw[i] = gamsor(k,q, TDMA(data, coeff))[0][i] / grid(x1, x2, data)
Ap[i] = Ae[i] + Aw[i]
# Define "b" array - including source term and boundary conditions
b = np.zeros((data(n, hw, he, TinfW, TinfE)[0], 1))
b.fill(gamsor(k,q, TDMA(data, coeff))[1][i])
b[0] = (gamsor(k,q, TDMA(data, coeff))[1][i]*(grid(x1,x2,data)/2) +
data(n, hw, he, TinfW, TinfE)[1]*
data(n, hw, he, TinfW, TinfE)[3])
b[-1] = (gamsor(k,q, TDMA(data, coeff))[1][i]*(grid(x1,x2,data)/2)+
data(n, hw, he, TinfW, TinfE)[2]*
data(n, hw, he, TinfW, TinfE)[4])
# Change values in first and last coefficients to reflect specified BC's
Aw[0] = 0
Ap[0] = Ae[0] + Aw[0] + data(n, hw, he, TinfW, TinfE)[1]
Ae[-1] = 0
Ap[-1] = Ae[0] + Aw[0] + data(n, hw, he, TinfW, TinfE)[2]
return Ae, Aw, Ap, b
#Step 5: SOLVE
def TDMA(data, coeff):
# Initialize "T" array - the solution array
T = np.zeros(data(n, hw, he, TinfW, TinfE)[0]);
# Initialize recursion functions as arrays
P = np.zeros(data(n, hw, he, TinfW, TinfE)[0]);
Q = np.zeros(data(n, hw, he, TinfW, TinfE)[0]);
## Step 1: evaluate at node 1
P[0] = (coeff(Ae, Aw, Ap, data, gamsor, grid)[0][0] /
coeff(Ae, Aw, Ap, data, gamsor, grid)[2][0])
Q[0] = (coeff(Ae, Aw, Ap, data, gamsor, grid)[3][0] /
coeff(Ae, Aw, Ap, data, gamsor, grid)[2][0])
## Step 2: sweep from node 2 to node (n-1) (python
## cell 1->(n-2) ) evaluating P and Q
for i in range(1, data(n, hw, he, TinfW, TinfE)[0]-1):
P[i] = ((coeff(Ae, Aw, Ap, data, gamsor, grid)[0][i]) /
(coeff(Ae, Aw, Ap, data, gamsor, grid)[2][i] -
coeff(Ae, Aw, Ap, data, gamsor, grid)[1][i]*P[i-1]))
Q[i] = ((coeff(Ae, Aw, Ap, data, gamsor, grid)[3][i] +
coeff(Ae, Aw, Ap, data, gamsor, grid)[1][i]*Q[i-1]) /
(coeff(Ae, Aw, Ap, data, gamsor, grid)[2][i] -
coeff(Ae, Aw, Ap, data, gamsor, grid)[1][i]*P[i-1]))
## Step 3: calculate for node n
P[data(n, hw, he, TinfW, TinfE)[0]-1] = 0
Q[data(n, hw, he, TinfW, TinfE)[0]-1] = (
(coeff(Ae, Aw, Ap, data, gamsor, grid)[3][data(n, hw, he,
TinfW, TinfE)[0]-1] + coeff(Ae, Aw, Ap, data, gamsor, grid)
[1][data(n, hw, he, TinfW, TinfE)[0]-1]*Q[data(n, hw, he, TinfW,
TinfE)[0]-2]) / (coeff(Ae, Aw, Ap, data, gamsor, grid)[2][data(n,
hw, he, TinfW, TinfE)[0]-1] - coeff(Ae, Aw, Ap, data, gamsor,
grid)[1][data(n, hw, he, TinfW, TinfE)[0]-1]*P[data(n, hw, he,
TinfW, TinfE)[0]-2]))
T[data(n, hw, he, TinfW, TinfE)[0]-1] = (
Q[data(n, hw, he, TinfW, TinfE)[0]-1] )
## Step 4: back fill, giving the temperatures
for i in range(data(n, hw, he, TinfW, TinfE)[0] - 2, -1, -1):
T[i] = P[i]*T[i+1] + Q[i]
return T
def Temps(TDMA):
Temps = copy.copy(TDMA(data,coeff)); Temps.fill(100)
return Temps
if __name__ == '__main__':
# Initialize variables
n = 0
hw = 0
he = 0
TinfW = 0
TinfE = 0
k = np.zeros(data(n, hw, he, TinfW, TinfE)[0])
q = np.zeros(data(n, hw, he, TinfW, TinfE)[0])
k0 = 0
beta = 20
T0 = 0
x1 = 0
x2 = 0
Ae = np.zeros(data(n, hw, he, TinfW, TinfE)[0])
Aw = np.zeros(data(n, hw, he, TinfW, TinfE)[0])
Ap = np.zeros(data(n, hw, he, TinfW, TinfE)[0])
b = 0
data(n, hw, he, TinfW, TinfE)
grid(x1, x2, data)
Temps(TDMA)
for _ in range(1):
gamsor(k, q)
coeff(Ae, Aw, Ap, data, gamsor, grid)
TDMA(data, coeff)
#print Temps(TDMA)
print "Nodal Temperature Distribution (C):", TDMA(data, coeff)
