Pass forward intermediate results during iterative optimization

To investigate a counter-current flow heat exchanger while considering temperature dependent physical properties (such as specific heat $$c_\textrm{p,i}$$, heat conductivity $$\lambda_\textrm{i}$$, viscosity $$\eta_\textrm{i}$$ ...) I am trying to solve a system of equations by finding its roots using scipy.optimize.fsolve().

The temperature dependance can be summarized as

$$\vartheta_\textrm{i,mean} = \frac{\vartheta_\textrm{i,in} + \vartheta_\textrm{i,out}}{2} \\ \eta_\textrm{i} = f(i, \vartheta_\textrm{i,mean}) \\ c_\textrm{p,i} = f(i, \vartheta_\textrm{i,mean}) \\ \lambda_\textrm{i} = f(i, \vartheta_\textrm{i,mean}) \\ \alpha_\textrm{i} = \frac{\textrm{Nu}_\textrm{i} \cdot \lambda_\textrm{i}}{d_\textrm{h}} \\$$

and is implemented with regression polynomials.

Based on these physical properties, dimensionless numbers are calculated (Reynold's number $$Re$$, Prandtl's number $$Pr$$ and Nusselt's number $$Nu$$:

$$\textrm{Re}_\textrm{i} = \frac{\dot{m}_\textrm{i} \cdot d_\textrm{h}}{A_\textrm{c} \cdot \eta_\textrm{i}} \\ \textrm{Pr}_\textrm{i} = \frac{\eta_\textrm{i} \cdot c_\textrm{p,i}}{\lambda_\textrm{i}} \\ \textrm{Nu}_\textrm{i} = \textrm{Pr}_\textrm{i}^{1/3} \cdot \textrm{Re}_\textrm{i}^{0.62} \\$$

As heat is transferred from a hot to a cold fluid stream, I need to evaluate the beforementioned equations for each fluids which are denoted by the general index $$i$$:

$$i = \{1, 2\}$$

The above equations are then combined with each other to calculate the thermal transmittance $$k$$:

$$k = \frac{1}{ \frac{1}{\alpha_\textrm{1} } + \frac{s_\textrm{w}}{\lambda_\textrm{w}} + \frac{1}{\alpha_\textrm{2}}}$$

The transferred heat $$\dot{Q}$$

$$\dot{Q} = \dot{m}_\textrm{1} \cdot c_\textrm{p,1} \cdot (\vartheta_\textrm{1,out} - \vartheta_\textrm{1,in}),$$

can then be used to define the system of equations to be finally solved:

$$0 = \vartheta_\textrm{2,mean} - \vartheta_\textrm{1,mean} - \Delta\vartheta_\textrm{mean} \\ 0 = \dot{m}_\textrm{2} \cdot c_\textrm{p,2} \cdot (\vartheta_\textrm{2,in} - \vartheta_\textrm{2,out}) - \dot{Q} \\ 0 = k \cdot A_\textrm{hex} \cdot \Delta\vartheta_\textrm{mean} - \dot{Q} \\$$

My Python implementation looks like this (with pseudo-code elements):

import numpy as np
from scipy.optimize import fsolve

# -- stream 1
m_flow_1 = 30
theta_1_in = 40.9
theta_1_out = 70    # initially guessed value for solver

# -- stream 2
m_flow_2 = 100
theta_2_in = 105.1
theta_2_out = 90    # initially guessed value for solver

# -- heat exchanger configuration
area_heatex = 2312.5e-6
lambda_wall = 15
wall_thickness = 0.75e-3

def alpha_1(*args):
# calculate alpha for fluid stream 1 here
pass

def alpha_2(*args):
# calculate alpha for fluid stream 2 here
pass

def heat_cap_1(theta):
# calculate heat capacity for fluid stream 1 here
pass

def heat_cap_2(theta):
# calculate heat capacity for fluid stream 2 here
pass

def nusselt_1(*args):
# calculate Nusselt's number for fluid stream 1 here
pass

def nusselt_2(*args):
# calculate Nusselt's number for fluid stream 2 here
pass

def thermal_transmittance(alpha_1, alpha_2, thickness, thermal_conductivity):
a = 1 / alpha_1
b = thickness / thermal_conductivity
c = 1 / alpha_2
return 1 / (a + b + c)

def theta_mean(theta_in, theta_out):
return (theta_in + theta_out) / 2

def delta_theta_mean(theta_2_in, theta_2_out, theta_1_in, theta_1_out):
return theta_mean(theta_2_in, theta_2_out) - theta_mean(theta_1_in, theta_1_out)

def func(z):
theta_1_out = z
theta_2_out = z

theta_1 = theta_mean(theta_1_in, theta_1_out)
nu_1 = nusselt_1(m_flow_1, theta_1)
alph_1 = alpha_1(nu_1, theta_1)

theta_2 = theta_mean(theta_2_in, theta_2_out)
nu_2 = nusselt_2(m_flow_2, theta_2)
alph_2 = alpha_2(nu_2, theta_2)

k = thermal_transmittance(alph_1, alph_2, wall_thickness, lambda_wall)
delta_t = delta_theta_mean(theta_2_in, theta_2_out, theta_1_in, theta_2_out)

cp_1 = heat_cap_1(theta_1)
cp_2 = heat_cap_2(theta_2)

q_dot = m_flow_1 * cp_1 * (theta_1_out - theta_1_in)

f = np.empty(2)
f = m_flow_2 * cp_2 * (theta_2_in - theta_2_out) - q_dot
f = k * area_heatex * delta_t - q_dot

return f

guessed = np.array([theta_1_out, theta_2_out])
solution = fsolve(func, guessed, full_output=True)
print(solution)

The function func containing the system of equations is passed to fsolve in order to find a valid solution/root.

As you can see from the above equations and the code implementation, new physical properties need to be calculated based on the temperature values calculated by the solver. However, this results in wrong physical properties as the values from each iteration step do not seem to represent the iteratively calculated temperatures but rather their relative deviation (?).

Is there a way to access the iteration results to calculate the required physical properties and dimensionless numbers inside the iteration function?