Some simple calculations of wort cooling with an internal coil

wort_cooling.py

import numpy as np
import matplotlib.pyplot as plt

from scipy.integrate import solve_ivp

# LMTD case
def f(t,y,params):

    p2, p3, Tin = params

    Tout = (1. - 1./np.exp(p2))*y + 1./np.exp(p2)*Tin
    return p3*(Tin - Tout)

# \infty-flow case
def ftin(t,y,params):

    p1, Tin = params
    return p1*(Tin - y[0])

# CSTR-case
def ftout(t,y,params):

    p1, p2, Tin = params

    return p1/(1.+p2)*(Tin - y[0])

# target event
def reaches_25(t,y):

        return y[0] - (273 + 25)

def main():

    phi = 0.01   # pretok [kg/s]
    phi = 200./3600.
    cpw = 4200. # specificna toplota vode [J/kg/K]
    mb = 15.    # masa piva [kg]
    cpb = 3849.   # specificna toplota piva

    U = 400. # toplotna prehodnost stene [W/m^2/K]
    A = 0.3  # povrsina spirale [m^2]

    Tin = 273.15 + 13 # vstopna temperature [K]

    tmax = 40
    t = (0,tmax*60) # time interval
    t_eval = np.linspace(t[0],t[1],100)

    y0 = [273.15 + 100] # zacetna temperatura piva [K]

    p1 = U*A/mb/cpb
    p2 = U*A/phi/cpw
    p3 = phi*cpw/mb/cpb

    flow = np.linspace(50./3600.,900./3600,60)
    t25 = []
    for flw in flow:
        p2flw = U*A/flw/cpw
        p3flw = flw*cpw/mb/cpb 
        sol = solve_ivp(fun=lambda t, y: f(t,y,(p2flw,p3flw,Tin)),t_span=t,y0=y0,t_eval=t_eval,events=reaches_25)
        t25.append(sol.t_events[0][0]/60.)

    plt.plot(flow*3600,t25,label="Time from 100 °C till 25 °C")
    plt.xlabel("Flow rate [kg/h]")
    plt.ylabel("Necessary time [min]")
    plt.legend()
    plt.show()


    sol = solve_ivp(fun=lambda t, y: f(t,y,(p2,p3,Tin)),t_span=t,y0=y0,t_eval=t_eval,events=reaches_25)
    sol_tin = solve_ivp(fun=lambda t, y: ftin(t,y,(p1,Tin)),t_span=t,y0=y0,t_eval=t_eval)
    sol_tout = solve_ivp(fun=lambda t, y: ftout(t,y,(p1,p2,Tin)),t_span=t,y0=y0,t_eval=t_eval)

    print("t = ", sol.t)
    print("t_eval = ", t_eval)
    print("y = ", sol.y)
    print(" reaches_25 ", sol.t_events[0][0]/60)


    tmin = t_eval/60
    
    plt.figure()
    plt.plot(tmin,sol_tin.y[0,:]-273,'-', color='C0', label=r"$\infty$-flow $T(t)$")
    
    Tout = (1-np.exp(-p2))*sol.y[0,:] + np.exp(-p2)*Tin
    plt.plot(tmin,sol.y[0,:]-273,'-',color='C1',label=r'LMTD $T(t)$')
    plt.plot(tmin,Tout-273,'--',color='C1',label=r'LMTD $T_{{out}}(t)$')

    
    Tout = (Tin + p2*sol_tout.y[0,:])/(1.+p2)
    plt.plot(tmin,sol_tout.y[0,:]-273,'-',color='C2',label=r"CSTR $T$")
    plt.plot(tmin,Tout-273,'--',color='C2',label=r'CSTR $T_{{out}}(t)$')



    plt.legend(loc=0)
    plt.hlines(25,tmin[0],tmin[-1],linestyle='dotted')
    plt.ylabel("Wort temperatue [°C]")
    plt.xlabel("Time [min]")
    plt.title(r"Wort chilling, $\phi = {}$ kg/h".format(phi*3600))
    plt.grid()
    plt.savefig("wort_chilling2.png")


    plt.show()

if __name__ == '__main__':
    main()