%Solution for Radial Heat Flow in a thermal pile using Laplace Transforms
%Nick Woodman, University of Southampton, June 2020

%Based on the model detailed in Report 2011:4, 'Radial Heat flow for a pipe in a
%borehole in ground using Laplace solutions Report on Mathematical
%background'. Johan Claesson, Department of Civil and
%Environmental Engineering, Chalmers University of Technology, Sweden.

%Please refer also to associated paper Loveridge, Woodman, Javed & Claesson (2022) A fast
%approximate method for simulating thermal pile heat exchangers. https://doi.org/10.1016/j.gete.2022.100368


function T = T_CJRM(Q,lamg,lamc,Cg,Cc,rb,rp,Cp,Rp,to,tseries);

%Calculate Thermal diffusivity
diffg=lamg/Cg; %Thermal diffusivity for the ground (m2/s)
diffc=lamc/Cc; %Thermal diffusivity for the pile (m2/s)

%Calculate thermal resistances
taup=rp/sqrt(diffc*to);
taub=rb/sqrt(diffc*to);
taug=rb/sqrt(diffg*to);


%Calculate kbt, kbp, kbb, kbg functions
kbt = (@(u) 4*lamc./( besselj(0,taup*u).*bessely(0,taub*u)-bessely(0,taup*u).*besselj(0,taub*u)));
kbp = (@(u) 4*lamc * (0.5*pi*taup.*u.*(besselj(1,taup*u).*bessely(0,taub*u)-bessely(1,taup*u).*besselj(0,taub*u))-1)./(besselj(0,taup*u).*bessely(0,taub*u)-bessely(0,taup*u).*besselj(0,taub*u)));
kbb = (@(u) 4*lamc * (0.5*pi*taub.*u.*(besselj(1,taub*u).*bessely(0,taup*u)-bessely(1,taub*u).*besselj(0,taup*u))-1)./(besselj(0,taup*u).*bessely(0,taub*u)-bessely(0,taup*u).*besselj(0,taub*u)));
kbg = (@(u) 2*pi * lamg* taug.*u.*(besselj(1,taug*u)-j*bessely(1,taug*u))./(besselj(0,taug*u)-j*bessely(0,taug*u)));

%Compute integral
Lu =(@(u) imag(-Q./ ((Cp*(-u.^2)/to)+ 1./(Rp+ 1./ (kbp(u)+1./((1./kbt(u))+1./(kbb(u)+kbg(u)))  ))  ) ));
count=0;

 for dum=tindexmin:indstep:tindexmax;
     t=tseries(dum)*to;
     count=count+1
    
%%Integrate with Gauss-Kronod quadrature
Tf=(2/pi)* quadgk(@(u) ((1-exp(-((u.^2)*t/to)))./u).*Lu(u), 0, inf);
T(count)=Tf;
 
 end

%concatenate output into columns
T=T';

end


%--Nomenclature-----------------------------------------------------------%
%Q Power per unit depth (W/m)
%lamg Thermal conducitivity of ground (W/mK)
%lamc Thermal conductivity of pile concrete grout(W/mK)
%lamp Thermal conductivity of pipe (W/mK)
%Cg Specific heat capacity of the ground (J/kgK)
%Cc Specific heat capacity of the pile (J/kgK)
%rb Radius of the pile. (m)
%rp Effective radius of the (single, central) water pipe (m)
%Cp Thermal capacity of the heat-transfer liquid within the pipe per unit depth (J/mK)
%Rp Steady-state thermal resistance for the pipe (mK/W)
%to Time constant, set to 3600 (s)
%tseries Simulation times (s)