Milk leaves the cow at about 37 °C and must be below 4 °C quickly, or bacteria double every half hour. Run the milking round: set the pump and coolant flow, choose well water or a glycol chiller, and see the plate-cooler effectiveness, the tank temperature, the bacterial count and the electricity bill.
Sweep one setting over a whole round and see which value meets the 4 °C / 2 h rule at the lowest running cost.
Compares the code with an independent numerical counterflow solver, the closed-form Newton cooling curve, a hand-computed reference and the cardinal-model growth law.
The science
Why cool fast. Bacterial growth rate depends strongly on temperature. The cardinal-temperature model (Rosso et al., 1993) gives the specific growth rate μ(T) = μopt·γ(T) between a minimum, optimum and maximum temperature:
γ(T) = (T−Tmax)(T−Tmin)² / { (Topt−Tmin) [ (Topt−Tmin)(T−Topt) − (Topt−Tmax)(Topt+Tmin−2T) ] }
d ln N/dt = μ(T), tdouble = ln 2 / μ
Two illustrative populations are tracked: mesophiles (Tmin 5, Topt 37, Tmax 46 °C, 0.6 h doubling at optimum) and psychrotrophs (−6, 27, 35 °C, 1.0 h). Parameters are typical textbook-style values, not a fit to a specific dairy; no lag phase is modelled (pessimistic) and no competition or spoilage thresholds are simulated.
Plate heat exchanger (effectiveness–NTU, counterflow).
Cm = ṁmcp,m, Cc = ṁccp,c, Cmin, Cr = Cmin/Cmax
NTU = UA / Cmin
ε = (1 − e−NTU(1−Cr)) / (1 − Cr e−NTU(1−Cr)) (Cr=1: ε = NTU/(1+NTU))
Q = ε Cmin (Tm,in − Tc,in), Tm,out = Tm,in − Q/Cm
The overall coefficient adds film resistances and fouling: 1/U = 1/hm + 1/hc + tw/k + Rf, with film coefficients scaling like (flow per channel)0.7. Area is (plates−1)×0.085 m². More plates give more area but also lower flow per channel, so returns diminish. Well water is cheap but needs about 2–3 litres per litre of milk; a glycol chiller reaches 0–2 °C but costs compressor electricity.
Bulk tank. A well-mixed tank, with new milk mixed in at the cooler outlet temperature:
M cp dT/dt = UAleak(Tamb − T) + Pagitator − Qcool
Newton cooling (no refrigeration): T(t) = Tamb + (T0 − Tamb) e−t/τ, τ = M cp/UAleak
The compressor removes Qcool = Pel·COP with COP = 0.42·Te/(Tc−Te) (about 42 % of the Carnot limit; Te = Tmilk−8 K, Tc = Tamb+12 K, in kelvin). A thermostat cuts in at set-point + 1 K and out at the set-point; the compressor only runs once milk covers the evaporator (60 kg).
Quality rule used for scoring. Rules differ by country; a common benchmark is milk at or below about 4 °C within 2 hours of the end of milking, never letting the blend rise above about 10 °C, and a plate count below about 100,000 CFU/mL for raw milk. The score and prices (electricity, 0.60 per m³ of well water) are teaching values, not advice.
Honest simplifications. Constant specific heats (milk 3.93, water 4.18, glycol mix 3.6 kJ/kg·K); steady-state exchanger at each time step (the pack's thermal mass is ignored); no pressure-drop model beyond a fixed pump power estimate; chiller treated as an ice-bank type with 12 kW capacity; jar milk stays at udder temperature; ideal mixing in the tank; no heat of agitation distribution; milk fat/viscosity effects ignored. Generic equipment, not any manufacturer's design.
Baked 3D scene. The room, equipment and cows are modelled in Blender 5.2 with procedurally baked albedo and normal maps and a baked floor occlusion map; Three.js r160 runs the live part.