p3_num_code

import numpy as np
import numpy.linalg as npla
import matplotlib.pyplot as plt


# geometry

r_f = 1.056/2
t_g = 0.008
t_c = 0.0864
r_c = 0.6224

# thermal properties

kf = 0.025
kg = 0.004
kc = 0.107

t_inlet = 269+273
t_outlet = 286+273
t_inf = (t_inlet+t_outlet)/2
t_inf_C = (269+286)/2

avg_heat = 50.3
max_heat = 111.5

avg_vol_heat = avg_heat*2/r_f
max_vol_heat = max_heat*2/r_f

h_inf_0 = 0.0484

# nodes

n_f = 10
n_g = 5
n_c = 5
n_total = n_f + n_g + n_c + 1

# spacing

df = r_f / n_f
dg = t_g / n_g
dc = t_c / n_c

# build list of radius at each point

rad_f = np.linspace(df, r_f, n_f)
rad_g = np.linspace(r_f + dg, t_g + r_f, n_g)
rad_c = np.linspace(rad_g[-1] + dc, r_c, n_c)

radii = np.append(rad_f, [rad_g, rad_c])
radii = np.append(np.array(0), radii)

# set id for each region

# 0 = fuel center node
# 1 = fuel inner node
# 2 = fuel/gap interface node
# 3 = gap inner node
# 4 = gap/clad interface node
# 5 = clad inner node
# 6 = clad/coolant interface node

region_id = np.zeros(n_total, dtype=int)

region_id[1:n_f] = 1
region_id[n_f] = 2
region_id[n_f + 1:n_f + n_g] = 3
region_id[n_f + n_g] = 4
region_id[n_f + n_g + 1:n_total - 1] = 5
region_id[n_total - 1] = 6

# set up functions for coefficients

# fuel nodes
# fuel inner backward (T_i-1)
f_back = lambda i: -(radii[i]-df/2)/radii[i]

# fuel inner forward (T_i+1)
f_forw = lambda i: -(radii[i]+df/2)/radii[i]

# fuel/gap interface nodes
# fuel/gap interface backward
fg_back = lambda i: kf/df*(radii[i]-df/2)

# fuel/gap interface central
fg_cent = lambda i: (-kf/df*(radii[i]-df/2)-kg/dg*(radii[i]+dg/2))

# fuel/gap interface forward
fg_forw = lambda i: (kg/dg*(radii[i]+dg/2))

# fuel/gap interface q
fg_q = lambda i: -(radii[i]*df-df**2/4)/2

# gap nodes
# gap inner backward
g_back = lambda i: (radii[i]-dg/2)/dg

# gap inner central
g_cent = lambda i: -2*radii[i]/dg

# gap inner forward
g_forw = lambda i: (radii[i]+dg/2)/dg

# gap/clad interface nodes
# gap/clad interface backward
gc_back = lambda i: kg/dg*(radii[i]-dg/2)

# gap/clad interface central
gc_cent = lambda i: (-kg/dg*(radii[i]-dg/2)-kc/dc*(radii[i]+dc/2))

# gap/clad interface forward
gc_forw = lambda i: kc/dc*(radii[i]+dc/2)

# clad inner nodes
# clad inner forward
c_forw = lambda i: (radii[i]+dc/2)/dc

# clad inner central
c_cent = lambda i: -2*radii[i]/dc

# clad inner backward
c_back = lambda i: (radii[i]-dc/2)/dc


# set up matrices

a_mat = np.zeros([n_total, n_total])
q_vec = np.zeros(n_total)

# fill the a matrix with coefficients
for i in range(n_total):
    # get the id of the current node
    node = region_id[i]
    # place coefficients according to the region id

    # node at center of fuel
    if node == 0:
        a_mat[i, i] = 1
        a_mat[i, i+1] = -1
        q_vec[i] = df**2/4/kf

    # node at inner fuel
    if node == 1:
        a_mat[i, i-1] = f_back(i)
        a_mat[i, i+1] = f_forw(i)
        a_mat[i, i] = 2
        q_vec[i] = df**2/kf

    # node at fuel/gap interface
    if node == 2:
        a_mat[i, i-1] = fg_back(i)
        a_mat[i, i+1] = fg_forw(i)
        a_mat[i, i] = fg_cent(i)
        q_vec[i] = fg_q(i)

    # node at inner gap
    if node == 3:
        a_mat[i, i-1] = g_back(i)
        a_mat[i, i+1] = g_forw(i)
        a_mat[i, i] = g_cent(i)

    # node at gap/clad interface
    if node == 4:
        a_mat[i, i-1] = gc_back(i)
        a_mat[i, i+1] = gc_forw(i)
        a_mat[i, i] = gc_cent(i)

    # node at inner clad
    if node == 5:
        a_mat[i, i-1] = c_back(i)
        a_mat[i, i+1] = c_forw(i)
        a_mat[i, i] = c_cent(i)

    # node at clad/coolant interface
    if node == 6:
        a_mat[i, i-1] = -kc/dc*(r_c-dc/2)
        a_mat[i, i] = h_inf_0*r_c+kc/dc*(r_c-dc/2)
        q_vec[i] = h_inf_0*r_c*t_inf_C


# solve the system
def solve_system(q_vector):
    temps = npla.solve(a_mat, q_vector)
    return temps


temp_profile = solve_system(q_vec)
fig, ax = plt.subplots()
ax.plot(radii, temp_profile)
plt.show()
fig.savefig("/users/1uw/documents/college/ne571/p3num.png")

plt.close("all")


def set_h(h_inf, q_vector):
    q = h_inf*r_c*t_inf_C
    q_vector[-1] = q
    return q_vector