In [1]:
%matplotlib inline

In [2]:
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sn

import pycollocation


Consider the following boundary value problem: \begin{align} \frac{d T_1}{d_A} =& -q,\ T_1(0) = T_{1,0} \\ \frac{d T_2}{d_A} =& -\frac{1}{2}q,\ T_2(A_{hx}) = T_{2,A_{hx}} \\ \end{align} where $q = (T_1 - T_2)U$.

This boundary value problem describes a counter current heat exchanger; a hot liquid enters a device and exchanges heat across a metal plate with a cold liquid traveling through the device in the opposite direction. Here, $T_1$ is the temperature of the hot stream, $T_2$ is the temperature of the cold stream, $q$ is the rate of heat transfer from the hot fluid to the cold fluid, $U$ is the overall heat transfer coefficient and $A_hx$ is the total area of the heat exchanger.

In [4]:
# provide numerical values for the parameters
params = {'T10': 130, 'T2Ahx': 70, 'U': 1.0}

def q(T1, T2, U):
return (T1 - T2) * U

def rhs(A, T1, T2, U, **params):
return [-q(T1, T2, U), -0.5 * q(T1, T2, U)]

def bcs_lower(A, T1, T2, T10, **params):
return [T1 - T10]

def bcs_upper(A, T1, T2, T2Ahx, **params):
return [T2 - T2Ahx]

In [5]:
problem = pycollocation.problems.TwoPointBVP(bcs_lower, bcs_upper, 1, 2, params, rhs)

In [8]:
def initial_mesh(domain, num, problem):
As = np.linspace(domain[0], domain[1], num)
T1s = np.repeat(0.5 * (problem.params['T10'] + problem.params['T2Ahx']), num)
return As, T1s, T1s

In [10]:
polynomial_basis = pycollocation.basis_functions.PolynomialBasis()
solver = pycollocation.solvers.Solver(polynomial_basis)

basis_kwargs = {'kind': 'Chebyshev', 'domain': [0, 5.0], 'degree': 15}

As, T1s, T2s = initial_mesh(basis_kwargs['domain'], 1000, problem)
T1_poly = polynomial_basis.fit(As, T1s, **basis_kwargs)
T2_poly = polynomial_basis.fit(As, T2s, **basis_kwargs)
initial_coefs = np.hstack([T1_poly.coef, T2_poly.coef])

solution = solver.solve(basis_kwargs, initial_coefs, problem)

In [11]:
T1_soln, T2_soln = solution.evaluate_solution(As)
plt.plot(As, T1_soln)
plt.plot(As, T2_soln)
plt.show()

In [12]:
T1_normalized_resids, T2_normalized_resids = solution.normalize_residuals(As)
plt.plot(As, np.abs(T1_normalized_resids))
plt.plot(As, np.abs(T2_normalized_resids))
plt.yscale('log')
plt.show()

In [ ]: