Modeling of terahertz heating of materials

From Traywiki
Jump to: navigation, search

Introduction

This software models the equilibrium temperature of a material when exposed to a THz (T-ray) beam. It is coded in Matlab and specifically analyses the case of a static sample of water. The parameters can be changed so that you can model any material. For further information, see: Optics Express, Vol. 18, No. 5, pp. 4727-4739, 2010.

Legal notice

This open source software is freely available. You may use it at will and adapt it. If you extend the software and wish us to put your new version here or have corrections to the existing version, please email us. Please note that your use of our software is entirely at your own risk and we do not offer any guaranteed support. It is available free and it is upto you to adapt it for your own purposes.

Citing an acknowlegment

If your use of this software results in a publication, our only request is that you acknowledge the source by citing the following reference:

[1] T. T. L. Kristensen, W. Withayachumnankul, P. U. Jepsen, and D. Abbott, "Modeling terahertz heating effects on water," Optics Express, Vol. 18, No. 5, pp. 4727-4739, 2010.

The code

function [Umax,conv,N,U] = thzheat(r,z,N,P,a,T0,nu,b,d,disP,bc,T)
% THZHEAT Calculates the temperature increase of a cylindrical sample of
% water. The model is valid in the frequency range 0 THz to 25 THz and in
% the temperature range 0 degC to 100 degC. The output Umax is in K at a
% particular location, (r,z), in the sample. Users can determine the
% temperature increase in an arbitrary dielectric sample by changing
% alpha and n, for the absorption coefficient and refractive index,
% in the function AbsCoef below.
%
% For further information see:
% T. T. Kristensen, W. Withayachumnankul, P. U. Jepsen and D. Abbott,
% "Modeling terahertz heating effects on water," Opt. Express 18, 4727-4739 (2010)
%
% Input parameters:
% r; radial coordinate, in metres (scalar or vector)
% z; z coordinate, in metres (scalar or vector)
% N; number of terms to include (scalar)
% P; input power, in Watts (scalar or vector)
% a; beam radius, in metres (scalar or vector)
% T0; initial temperature, in Kelvin (scalar or vector)
% nu; frequency, in Hz (scalar or vector)
% b; sample radius, in metres (scalar or vector)
% d; sample thickness, in metres (scalar or vector)
%
% disP (OPTIONAL); dissipated power density,
% 1={constant}, 2={exponential} (DEFAULT)
% bc (OPTIONAL); boundary condition for the side of the disk
% 1={bc b.1}, 2={bc b.2}, 3={the term C_n is neglected} (DEFAULT)
% T (OPTIONAL); transmittance at the interface
% 0={exclude transmittance} (DEFAULT), 1={include transmittance}
% Output parameters:
% Umax; temperature increase, in K
% conv (OPTIONAL); accuracy of the value of the temperature increase
% N (OPTIONAL); number of terms used
% U (OPTIONAL); vector containing the value for each term in the sum
%
% Example:
% r = [0:0.01:1].*1e-3;
% [Umax,conv,N,U]=thzheat(r,0,1000,1e-3,0.25e-3,273.15+25,1e12,50e-3,15e-3,2,3,0);
% plot(r, Umax); xlabel('Radius [m]'); ylabel('Temperature increase [K]');
%
% Author: Torben T. L. Kristensen
% Date: January 11, 2010
if nargin<12
T=0;
end
if nargin<11
bc=3;
end
if nargin<10
disP=2;
end
%
% The number of terms need to obey the following equation, if too many
% terms are chosen N is automatically decreased to the maximum value.
%
maxN=floor((697*d./a*2/pi+1)/2);
if N>maxN
N=max(maxN);
disp(['N too large, new N = ',num2str(N)]);
end
%
% If any of the variables are a vector, a meshgrid is created between
% the variable and n.
n=1:N;
if length(r)>1
[r,n]=meshgrid(r,n);
elseif length(z)>1
[z,n]=meshgrid(z,n);
elseif length(P)>1
[P,n]=meshgrid(P,n);
elseif length(a)>1
[a,n]=meshgrid(a,n);
elseif length(T0)>1
[T0,n]=meshgrid(T0,n);
elseif length(nu)>1
[nu,n]=meshgrid(nu,n);
elseif length(b)>1
[b,n]=meshgrid(b,n);
elseif length(d)>1
[d,n]=meshgrid(d,n);
end
%
% The thermal conductivity, the absorption coefficient and the refractive
% index is calculated using the two subfunctions written below
%
k0=TermCon(T0);
[alpha,index]=AbsCoef(nu,T0);
c0=2.99792458e8;
index_complex=index+1i*alpha*c0./(4*pi*nu);
if ~T
index_complex=1;
end
% The temperature increase is calculated
mm=disP-1;
m2=zeros(size(n));
m2idx=find(r>a);
m2(m2idx)=1;
m3=abs(m2-1);
G=(mm*alpha+abs(mm-1)*1./d)./(pi*a.^2.*k0);
pn=(2*n-1)*pi/2./d;
pna=pn.*a;
pnb=pn.*b;
pnr=pn.*r;
gn=2*G./d.*(mm*alpha-(-1).^(n).*pn.*exp(-mm*alpha.*d))./((mm*alpha).^2+pn.^2);
An=1./(besseli(1,pna).*besselk(0,pna)+besseli(0,pna).*besselk(1,pna));
Bn=(-1).^m2.*besselk(m3,pna.*m3+pnr.*m2).*besseli(m2,pna.*m2+pnr.*m3);
if bc==3
U1=cos(pn.*z).*gn./pn.^2.*(m3-An.*Bn);
else
m=bc-1;
C2stop=floor((697*d./b*2/pi+1)/2);
Cn=zeros(size(n));
C2idx=find(n<C2stop);
Cn(C2idx)=(-1).^(m).*besseli(1,pna(C2idx)).*besselk(m,pnb(C2idx))./besseli(m,pnb(C2idx)).*besseli(0,pnr(C2idx).*m2(C2idx));
U1=cos(pn.*z).*gn./pn.^2.*(m3-An.*(Bn+Cn));
end
U1=U1.*P;
if size(U1,1)==1, U1=U1'; end
conv=U1(N,:);
tau=abs(4*index_complex./(1+index_complex).^2);
U=tau*cumsum(U1,1);
Umax=max(U);
end
% Function for the thermal conductivity
function k0=TermCon(T0)
k0=(-1.48445+4.12292/298.15*T0-1.63866/298.152^2*T0.^2)*0.6065;
% Costumized value for the thermal conductivity can be typed in here:
% k0=xxx
end
% Function for the absorption coefficient and the refractive index
function [alpha,index]=AbsCoef(nu,T0)
T=T0-273.15;
c0=2.99792458e8;
% Parameters for the three relaxations terms
a=[79.23882 3.815866 1.634967];
b=[0.004300598 0.01117295 0.006841548];
c=[1.382264E-13 3.510354E-16 6.30035E-15];
d=[652.7648 1249.533 405.5169];
Tc=133.1383;
% Resonance parameters
p=[0.8379692 -0.006118594 -0.000012936798 4235901000000 -14260880000 273815700 -1246943 9.618642E-14 1.795786E-16 -9.310017E-18 1.655473E-19...
0.6165532 0.007238532 -0.00009523366 15983170000000 -74413570000 497448000 2.882476E-14 -3.142118E-16 3.528051E-18];
es=87.9144-0.404399*T+9.58726E-4*T.^2-1.32802e-6*T.^3;
Delta{1}=a(1)*exp(-b(1)*T);
Delta{2}=a(2)*exp(-b(2)*T);
Delta{3}=a(3)*exp(-b(3)*T);
tau{1}=c(1).*exp(d(1)./(T+Tc));
tau{2}=c(2).*exp(d(2)./(T+Tc));
tau{3}=c(3).*exp(d(3)./(T+Tc));
Delta{4}=p(1)+p(2)*T+p(3)*T.^2;
f{1}=p(4)+p(5)*T+p(6)*T.^2+p(7)*T.^3;
tau{4}=p(8)+p(9)*T+p(10)*T.^2+p(11)*T.^3;
Delta{5}=p(12)+p(13)*T+p(14)*T.^2;
f{2}=p(15)+p(16)*T+p(17)*T.^2;
tau{5}=p(18)+p(19)*T+p(20)*T.^2;
K1={f{1}+nu f{1}-nu f{2}+nu f{2}-nu};
C1{1}=tau{1}.^2.*Delta{1}./(1+(2*pi*nu.*tau{1}).^2);
C1{2}=tau{2}.^2.*Delta{2}./(1+(2*pi*nu.*tau{2}).^2);
C1{3}=tau{3}.^2.*Delta{3}./(1+(2*pi*nu.*tau{3}).^2);
C2{1}=nu.*K1{1}./(1+(2*pi*tau{4}.*K1{1}).^2);
C2{2}=nu.*K1{2}./(1+(2*pi*tau{4}.*K1{2}).^2);
C2{3}=nu.*K1{3}./(1+(2*pi*tau{5}.*K1{3}).^2);
C2{4}=nu.*K1{4}./(1+(2*pi*tau{5}.*K1{4}).^2);
C3{1}=(2*pi*tau{4}).^2.*Delta{4}/2;
C3{2}=(2*pi*tau{5}).^2.*Delta{5}/2;
C4{1}=tau{1}.*Delta{1}./(1+(2*pi*nu.*tau{1}).^2);
C4{2}=tau{2}.*Delta{2}./(1+(2*pi*nu.*tau{2}).^2);
C4{3}=tau{3}.*Delta{3}./(1+(2*pi*nu.*tau{3}).^2);
C5{1}=1./(1+(2*pi*tau{4}.*K1{1}).^2);
C5{2}=1./(1+(2*pi*tau{4}.*K1{2}).^2);
C5{3}=1./(1+(2*pi*tau{5}.*K1{3}).^2);
C5{4}=1./(1+(2*pi*tau{5}.*K1{4}).^2);
C6{1}=pi*nu.*tau{4}.*Delta{4};
C6{2}=pi*nu.*tau{5}.*Delta{5};
C0=2*pi*nu;
e1=es-C0.^2.*(C1{1}+C1{2}+C1{3})-C3{1}.*(C2{1}-C2{2})-C3{2}.*(C2{3}-C2{4});
e2=C0.*(C4{1}+C4{2}+C4{3})+C6{1}.*(C5{1}+C5{2})+C6{2}.*(C5{3}+C5{4});
% Calculation of the absorption coefficient.
alpha=4*pi*nu./c0.*sqrt((e1.^2+e2.^2).^(1/2)/2-e1/2);
% Calculation of the refractive index.
index=sqrt((e1+sqrt(e1.^2+e2.^2))/2);
% Customized value for the absorption coefficient [1/cm] and refractive index can
% be typed in here:
% alpha=xxx
% index=xxx
end

See Also

External Links

Back

Retrieved from "https://www.eleceng.adelaide.edu.au/groups/thz/wiki/index.php?title=Modeling_of_terahertz_heating_of_materials&oldid=3096"