c c hotspot - model heat dissipation by conduction and radiation c in two dimensions under a localized heat load c c Richard T. Jones, July 23 2000 c c Modification history: c * August 17, 2006 -Richard Jones, G. Yang c - Introduced thermal relief due to conduction down the mounting wires c to the model. This drain is applied to the four corners of the c crystal and depends on the local temperature in the corners. c * October 16, 2004 -Richard Jones c - Added emissivity to the model, before I was assuming that diamond is c a blackbody which is a very bad assumption! c - Put in the physical heat capacity and changed to a physical time c step so that the simulation can cover the case of a pulsed beam. c c-------------------------------- c Model c c The model consists of a thin film characterized by c 1) the film thickness z (cm) c 2) bulk density rho of diamond (g/cm^3) c 3) heat capacity Cp (Joules/K/g) c 4) thermal conductivity kappa (Watts/K/cm) c 5) thermal emissivity ems of diamond (dimensionless) c 6) the ambient temperature T_0 of the surrounding space (Kelvin) c 7) Steffan-Boltzmann radiation constant SBc (Watts/cm^2/K^4) c 8) diameter and length of four supporting wires c 9) thermal conductivity of the supporting wires (Watts/cm*K) c The driving term in the model is the heating function h(x,y) that is c expressed in Watts/area. The basic unit for length is the cm, sorry c for mixing cgs and mks units. Time in this model is measured in c seconds, current in Amperes, energy in Joules and power in Watts, and c temperatures in degrees Kelvin. The equation of motion is as follows. c c dT 1 4 4 _2 c -- = -------- ( h(x,y) - 2 ems SBc (T - T ) + z kappa V T ) c dt z rho Cp 0 c c The funny V with a _2 overhead is the LaPlacian. Different choices c for the boundary conditions are provided by the diffuseX() functions. c The heat drain due to conductivity of the support wires is considered c to be a part of the boundary conditions. c Procedure c c 1) call initialize() to initialize the thermal grid. This will set c the temperature everywhere to a uniform value given by T_0. c 2) call hotspot(current) to set the total power level of the heater. c This can be called at any time to change the heating term. c 3) call diffuse(nsteps,dt) to advance nsteps in time steps of size dt. c The state of the thermal grid can be sampled at any time. subroutine hotspot(pulse_energy, pulse_duration) implicit none real pulse_energy, pulse_duration include 'hotspot.inc' data thickness,width,height /20.e-4,0.6,0.6/ data beamsizeX,beamsizeY,hpower/0.0830,0.0370,25.e-3/ real BetheBloch external BetheBloch real area area = 2*pi_*beamsizeX*beamsizeY hpower = pulse_energy/pulse_duration print *, 'Power on crystal in Watts is',hpower hpower = hpower/area end subroutine initialize() implicit none include 'hotspot.inc' integer ix,iy vector Tgrid(30,30) vector time(1) time(1) = 0 do ix=1,NX do iy=1,NY Tgrid(ix,iy) = Tambient enddo enddo end subroutine initialize1(Tinitial) implicit none real Tinitial include 'hotspot.inc' integer ix,iy vector Tgrid(30,30) vector time(1) time(1) = 0 do ix=1,NX do iy=1,NY Tgrid(ix,iy) = Tinitial enddo enddo end subroutine diffuse(nsteps,dtstep) implicit none integer nsteps real dtstep include 'hotspot.inc' vector Tgrid(30,30) vector time(1) integer i,ix,iy do ix=1,NX do iy=1,NY Tgrid_(ix,iy) = Tgrid(ix,iy) enddo enddo time_ = time(1) do i=1,nsteps call diffuse1(dtstep) enddo time(1) = time_ do ix=1,NX do iy=1,NY Tgrid(ix,iy) = Tgrid_(ix,iy) enddo enddo end subroutine diffuse1(dtstep) ! with thermal clamps at periphery implicit none real dtstep include 'hotspot.inc' real powerdensity external powerdensity double precision T(NX,NY),dT real x,y integer ix,iy do ix=1,NX do iy=1,NY T(ix,iy) = Tgrid_(ix,iy) enddo enddo time_ = time_+dtstep do ix=2,NX-1 x = ((ix-0.5)/NX-0.5)*width do iy=2,NY-1 y = ((iy-0.5)/NY-0.5)*height dT = 1./(rhoCp*thickness) * ( powerdensity(x,y) + - 2*ems*SBc*(T(ix,iy)**4 - Tambient**4) + + kappa*thickness * ( + (T(ix-1,iy) + T(ix+1,iy) - 2*T(ix,iy)) + /(width/NX)**2 + + (T(ix,iy-1) + T(ix,iy+1) - 2*T(ix,iy)) + /(height/NY)**2 + ) + ) Tgrid_(ix,iy) = Tgrid_(ix,iy) + dT*dtstep enddo enddo end subroutine diffuse2(dtstep) ! insulated, only radiative cooling implicit none real dtstep include 'hotspot.inc' real powerdensity external powerdensity double precision T(NX,NY),dT real x,y integer ix,iy vector dTmax(1) do ix=1,NX do iy=1,NY T(ix,iy) = Tgrid_(ix,iy) enddo enddo time_ = time_+dtstep dTmax(1) = 0 do ix=2,NX-1 x = ((ix-0.5)/NX-0.5)*width do iy=2,NY-1 y = ((iy-0.5)/NY-0.5)*height dT = 1./(rhoCp*thickness) * ( powerdensity(x,y) + - 2*ems*SBc*(T(ix,iy)**4 - Tambient**4) + + kappa*thickness * ( + (T(ix-1,iy) + T(ix+1,iy) - 2*T(ix,iy)) + /(width/NX)**2 + + (T(ix,iy-1) + T(ix,iy+1) - 2*T(ix,iy)) + /(height/NY)**2 + ) + ) Tgrid_(ix,iy) = Tgrid_(ix,iy) + dT*dtstep dTmax(1) = max(dTmax(1),dT*dtstep) enddo Tgrid_(ix,1) = Tgrid_(ix,2) Tgrid_(ix,NY) = Tgrid_(ix,NY-1) enddo do iy=1,NY Tgrid_(1,iy) = Tgrid_(2,iy) Tgrid_(NX,iy) = Tgrid_(NX-1,iy) enddo end subroutine diffuse3(dtstep) ! insulated along vertical boundaries implicit none real dtstep include 'hotspot.inc' real powerdensity external powerdensity double precision T(NX,NY),dT real x,y integer ix,iy do ix=1,NX do iy=1,NY T(ix,iy) = Tgrid_(ix,iy) enddo enddo time_ = time_+dtstep do ix=2,NX-1 x = ((ix-0.5)/NX-0.5)*width do iy=2,NY-1 y = ((iy-0.5)/NY-0.5)*height dT = 1./(rhoCp*thickness) * ( powerdensity(x,y) + - 2*ems*SBc*(T(ix,iy)**4 - Tambient**4) + + kappa*thickness * ( + (T(ix-1,iy) + T(ix+1,iy) - 2*T(ix,iy)) + /(width/NX)**2 + + (T(ix,iy-1) + T(ix,iy+1) - 2*T(ix,iy)) + /(height/NY)**2 + ) + ) Tgrid_(ix,iy) = Tgrid_(ix,iy) + dT*dtstep enddo enddo do iy=1,NY Tgrid_(1,iy) = Tgrid_(2,iy) Tgrid_(NX,iy) = Tgrid_(NX-1,iy) enddo end subroutine diffuse4(dtstep) ! metal wires connected at four corners implicit none real dtstep include 'hotspot.inc' real powerdensity external powerdensity double precision T(NX,NY),dT double precision dpower real x,y integer ix,iy do ix=1,NX do iy=1,NY T(ix,iy) = Tgrid_(ix,iy) enddo enddo time_ = time_+dtstep do ix=2,NX-1 x = ((ix-0.5)/NX-0.5)*width do iy=2,NY-1 y = ((iy-0.5)/NY-0.5)*height dT = 1./(rhoCp*thickness) * ( powerdensity(x,y) + - 2*ems*SBc*(T(ix,iy)**4 - Tambient**4) + + kappa*thickness * ( + (T(ix-1,iy) + T(ix+1,iy) - 2*T(ix,iy)) + /(width/NX)**2 + + (T(ix,iy-1) + T(ix,iy+1) - 2*T(ix,iy)) + /(height/NY)**2 + ) + ) Tgrid_(ix,iy) = Tgrid_(ix,iy) + dT*dtstep enddo Tgrid_(ix,1) = Tgrid_(ix,2) Tgrid_(ix,NY) = Tgrid_(ix,NY-1) enddo do iy=1,NY Tgrid_(1,iy) = Tgrid_(2,iy) Tgrid_(NX,iy) = Tgrid_(NX-1,iy) enddo do ix=2,NX-1,NX-3 do iy=2,NY-1,NY-3 dpower = (Tgrid_(ix,iy) - Tambient)/lengthWire + *kappaWire*(pi_*(diameterWire**2)/4) c dT = dpower * (height/NY) / (thickness * width * kappa) dT = dpower * NX*NY / (height * width * thickness * rhoCp) Tgrid_(ix,iy) = Tgrid_(ix,iy) - dT*dtstep enddo enddo end real function BetheBloch(E) ! in eV implicit none real E include 'hotspot.inc' real Tmax,Tcut,beta,gamma,del real hbwp Tmax = E/2 Tcut = min(Tmax,15.e6*density*thickness) gamma = 1+E/me_c2 beta = sqrt(1-1/gamma**2) hbwp = 28.816*sqrt(density*Z_/A_) del = 2*(log(hbwp/I_)+log(beta*gamma)-0.5) BetheBloch = K_*(Z_/A_)/beta**2 * ( + 0.5*log(2*me_c2*Tcut*(beta*gamma/I_)**2) + - (beta**2)/2*(1+Tcut/Tmax) + - del/2 + ) end real function powerdensity(x,y) ! in W/cm^2 implicit none real x,y include 'hotspot.inc' powerdensity = hpower * exp(-0.5*( + (x/beamsizeX)**2 + (y/beamsizeY)**2)) end