I am trying to solve coupled heat conduction and mass diffusion for a spherical particle.
I am solving the PDE’s using Finite difference method , with 4th accuracy for first and second order equations.
I am discretizing the spherical particle in 20 points .
num_points = 20
delta = radius/num_points.
I am creating an array with 20 elements as initial state for solving the equations.
X0 = np.zeros((num_points,))
I am trying the to solve the following equations:
1. For diffusion term.
def points_with_conc(T, dt, t)
if i != num_points-1: # For nodes 1 to 19 , from centre to second last node.
diff1 = Deff*(16*T[i+1]-T[i+2]-30*T[i]+16*T[i-1]-T[i-2])/12*delta**2
diff2 = Deff*(b-1)*(8*T[i+1]-T[i+2]-8*T[i-1]+T[i-2])/(12/delta/r_act)
diff = (diff1+diff2)/eps # Effective diffusion term
dXdt = get_dXdt(X, T[i], T[i+num_points]) # Calling the mass source function
reac = (dXdt*rho_s/gamma/44e-3/eps) # Effective Reaction rate term
else: # Node-20 , i.e At the Surface
cs = solve_cs(T[i], cg, eps) # Calling the surface
diff = (Deff*((16*T[i+1]-cs-30*T[i]+16*T[i-1]-T[i-2])/12*delta**2 + (b-1)*(8*T[i+1]-cs-8*T[i-1]+T[i-2])/(12/delta/r_act)/eps))
dXdt = get_dXdt(X, T[i], T[i+num_points])
reac = (dXdt*rho_s/gamma/44e-3/eps)
2. For Heat conduction term
if i != num_points-1: # For nodes 1 to 19 , from centre to second last node.
#Heat Diffusion Term - (∂T)/(∂t)=(α*(∂^2 T)/(∂r^2 ) + 2/r α*(∂T)/(∂r) + ∂α/∂r*(∂T)/(∂r))+Q
ai = lambda_i/rho_i/csolid # Effective Thermal Diffusivity
diff1 = ai*(16*T[i+num_points+1]-T[i+num_points+2]-30*T[i+num_points]+16*T[i+num_points-1]-T[i+num_points-2])/12*delta**2 # diffusion terms as shown in the HD equation above
diff2 = ai*(b-1)*(8*T[i+num_points+1]-T[i+num_points+2]-8*T[i+num_points-1]+T[i+num_points-2])/12/delta/r_act
diff = diff1+diff2 # Average diffusion term
dXdt = get_dXdt(X, T[i], T[i+num_points]) # Calling the mass source function
reac = dXdt*delta_rh/csolid/gamma*rho_s/rho_i # Effective reaction rate term
else: # Node-20 , i.e At the Surface
rho_i = X*rho_carb+(1-X)*rho_calc
lambda_i = (lambda_calc*lambda_carb)/(V_carb*lambda_carb + V_calc*lambda_calc)
ai = lambda_i/rho_i/csolid
Ts = solve_Ts(T[i+num_points], Tg, lambda_i) # Calling of surface temp function
diff1 = ai*(16*T[i+num_points+1]-Ts-30*T[i+num_points]+16*T[i+num_points-1]-T[i+num_points-2])/12*delta**2 # diffusion terms as shown in the HD equation above
diff2 = ai*(b-1)*(8*T[i+num_points+1]-Ts-8*T[i+num_points-1]+T[i+num_points-2])/12/delta/r_act
diff = float(diff1+diff2)
dXdt = get_dXdt(X, T[i], T[i+num_points])
reac = float(dXdt*delta_rh/csolid/gamma*rho_s/rho_i)
I am getting the following error, when i execute this code.
diff1 = ai*(16*T[i+num_points+1]-T[i+num_points+2]-30*T[i+num_points]+16*T[i+num_points-1]-T[i+num_points-2])/12*delta**2 # diffusion terms as shown in the HD equation above
IndexError: index 40 is out of bounds for axis 0 with size 40
Need some help solving this.
Thanks in advance.