Transient Thermal Lensing and Thermoelastic Abberations of Mirrors

The formulation is limited to right circular cylindrical mirrors with homogeneous surface or in-depth absorption:
Patrice Hello and Jean-Yves Vinet, "Analytical models of thermal aberrations in massive mirrors heated by high power laser beams", J. Phys. France, 51, (1990), 1267-1282
"Analytical models of transient thermoelastic deformations of mirrors heated by high power cw laser beams", J. Phys. France, 51, (1990), 2243-2261.

Dennis Coyne

Version History

Version 1.0 Initial

Version 2.0 restructured to have all inputs at beginning of the notebook and made transient calculations inactive for now.

Version 2.1 (a) changed the Zernike expansion to have the optical zone radius, b, be an input variable so that the Zernikes can be calculated for a nominal optical zone  radius b  and for the optic radius, (b) calculated the errors associated with the Zernike expansion

Version 2.2  Corrected "t" references to "tau" in the transcendental solutions for ur and vr. Implemented Zernike fit comparison to calculated surface/wavefront more uniformly in the notebooks (was only used in some and in only parts of the notebook.) Discovered and corrected error in the transcendental solutions ur and vr as used in the "analytic_ITM" notebook. Moved the calculation of sp and cp to after calculation of ur and vr. Corrected the calculation of the OPD at the center from the Zernike decomposition (changed from a sum of absolute amplitudes to an alternating signed summation).

Version 2.3 Applied to LIGO Optical Contamination Facility (OCF). See T970212

Preliminaries & Setup

Initialization

Off[General::"spell"]
Off[General::"spell1"]

Input Parameters

Optical Contamination Cavity Optic

a = optic radius, m
h = optic thickness, m
k = optic thermal conductivity, W/m/K
w = Gaussian beam waist, m
P = incident power, W
Tx = ambient temperature, K

eps = coating absorptivity
alpha = in-depth absorption, 1/m
Note: If the coating absorptivity == 1 or the in-depth absorption == 1/h , then the incident power == absorbed power.

ro = mass density, kg/m^3
c = specific heat, J/kg/K
dndT = change in refraction index with temperature, 1/K
b = radius(?) for Zernike decomposition, m
lambda = first Lame coefficient, J/m^3
mu = second Lame coefficient, J/m^3

a=0.0254/2;
h=0.0064;
k=1.38;
w=0.0004;
P=1;
Tx=293;

required for surface heating:

eps=1;

required for in-depth heating:

alpha=(1/h);

required for transient analysis:

ro=2202;
c=745;

required for thermal lensing:

dndT=(11.8/(10^6));

required for Zernike decomposition:

b=0.05;

required for thermoelastic analysis:

lambda=15.6 (10^9);
mu=31.3 (10^9);

ni=5.91 (10^4);

Constants

s=UnitConvert[Quantity["StefanBoltzmannConstant"],"Watts"/"Kelvins"^4/"Meters"^2][[1]]

Derived Parameters

tau = reduced radiation constant
tc= thermal time constant (for the entire optic), sec

tau=((4 s (Tx^3) a)/k);

tc=((ro c (a^2))/k)

(tc/3600)

Axisymmetric Zernikes

Z2q[q_,r_]:=(∑_(i=0)^q)(((((-1))^i) (2 q-i)! (r^(2 q-2 i)))/(i! (q-i)! (q-i)!));

Surface Absorption

Steady-State Thermal Field Solution

1st Transcendental Equation Solution

Equation Plot

Finds roots of the equation u J1(u) = c J0(u).  The user defines the constant, c, the root search interval, du, and the upper limit of u, upmax.

du=0.05;
upmax=200;
c=tau;

n=(upmax/du)+1
up=Table[i du,{i,n}];

e=N[up BesselJ[1,up]-c BesselJ[0,up]];
ListPlot[e]

Solution Method:

If the parameters are unchanged, then one can skip the calculation of the transcendental equation solution and go directly to the list of roots below.

es=Rest[e RotateRight[e]];
esw=Select[es,#1<0&];
npos=Dimensions[esw][[1]];
pos=Table[0,{npos}];
Do[pos[[i]]=Position[es,esw[[i]]],{i,1,npos}];
pos=Flatten[pos];
nu=Dimensions[pos][[1]];

ur=Table[0,{nu}];
Do[umin=up[[pos[[i]]]];umax=up[[pos[[i]]+1]];ustart=((umax+umin)/2);ur[[i]]=u/.FindRoot[u BesselJ[1,u]==c BesselJ[0,u],{u,ustart,umin,umax},AccuracyGoal->7,PrecisionGoal->7],{i,1,nu}];

ur

N[(ur-RotateRight[ur])/Pi]

mmax=Dimensions[ur][[1]];
zetam=ur;

Roots:

If the parameters are unchanged, then skip calculation above and evaluate this list of roots.

mmax=;
zetam={};

Coefficients of the Dini Series Representation of the Incident Radiation

zetam2=(zetam^2);
gamma=((zetam h)/(2 a));

pm=(P zetam2 Exp[-(((w^2) zetam2)/(8 (a^2)))])/(π (a^2) ((zetam2+(tau^2)) ((BesselJ[0,zetam])^2)));

Steady-State Temperature Distribution for Radiant Surface Heating

T[r_,z_]:=(∑_(m=1)^mmax)((pm[[m]] a) Exp[-((h zetam[[m]])/(2 a))] ((zetam[[m]]-tau) Exp[-((zetam[[m]] (h-z))/a)]+(zetam[[m]]+tau) Exp[-((zetam[[m]] z)/a)]) BesselJ[0,((zetam[[m]] r)/a)])/(k ((((zetam[[m]]+tau))^2)-(((zetam[[m]]-tau))^2) Exp[-((2 zetam[[m]] h)/a)]));

temp=Table[T[i,j],{i,0,a,(a/50)},{j,-(h/2),(h/2),(h/25)}];

Max[temp]

Min[temp]

ListPlot3D[temp,PlotRange->All,MeshRange->{{-(h/2),(h/2)},{0,a}},ViewPoint->{3,1,.5}]

ListContourPlot[Reverse[Transpose[temp]],PlotRange->All,Contours->10,AspectRatio->h/a,PlotLegends->Automatic]

Compare to ANSYS nonlinear, steady - state thermal reponse for the OCF mirror with 1 Watt absorbed power :
    max temperature = 813 C

Transient Thermal Field Solution

2nd Transcendental Equation Solution

Solution Method:

Finds roots of the equation u = t  Cot[u h/(2 a)].  The user defines the constant, t, the root search interval, du, and the upper limit of u, upmax.

dui=0.0001;
duj=((2 a π)/h);
upmax=120;
offset=0;

nj=Floor[(upmax/duj)+1];
ni=Floor[(duj/(10 dui))+1];
up=Flatten[Table[i dui+j duj+offset,{j,0,nj},{i,-ni,ni}]];

Dimensions[up]

e=N[up-tau Cot[((up h)/(2 a))]];
n=Dimensions[e][[1]];

ListPlot[Transpose[Partition[Join[up,e],n]]]

es=Rest[e RotateRight[e]];
esw=Select[es,#1<0&];
esw=Transpose[Partition[esw,2]][[1]];
npos=Dimensions[esw][[1]];
pos=Table[0,{npos}];
Do[pos[[i]]=Position[es,esw[[i]]],{i,1,npos}];
pos=Flatten[pos];
nu=Dimensions[pos][[1]];

pos

npos

ur=Table[0,{nu}];
Do[umin=up[[pos[[i]]]];umax=up[[pos[[i]]+1]];ustart=((umax+umin)/2);ur[[i]]=u/.FindRoot[u==tau Cot[((u h)/(2 a))],{u,ustart,umin,umax},MaxIterations->100,DampingFactor->0.8],{i,1,nu}];

ur

N[((ur-RotateRight[ur])/duj)]

ur2=ur^2;

Roots:

imax=;
ur={};
ur2=(ur^2);

3rd Transcendental Equation Solution

Finds roots of the equation u = -t  Tan[v h/(2 a)].  The user defines the constant, t, the root search interval, du, and the upper limit of v, vpmax.

Equation Plot:

Solution Method:

dvi=0.002;
dvj=((2 a π)/h);
vpmax=120;
offset=(dvj/2);

nj=Floor[(vpmax/dvj)+1];
ni=Floor[(dvj/(10 dvi))+1];
vp=Flatten[Table[i dvi+j dvj+offset,{j,0,nj},{i,-ni,ni}]];

Dimensions[vp]

e=N[vp+tau Tan[((vp h)/(2 a))]];
n=Dimensions[e][[1]];

ListPlot[e]

es=Rest[e RotateRight[e]];
esw=Select[es,#1<0&];
esw=Transpose[Partition[Drop[esw,1],2]][[1]];
npos=Dimensions[esw][[1]];
pos=Table[0,{npos}];
Do[pos[[i]]=Position[es,esw[[i]]],{i,1,npos}];
pos=Flatten[pos];
nu=Dimensions[pos][[1]];

npos

vr=Table[0,{nu}];
Do[vmin=vp[[pos[[i]]]];vmax=vp[[pos[[i]]+1]];vstart=((vmax+vmin)/2);vr[[i]]=v/.FindRoot[v==-tau Tan[((v h)/(2 a))],{v,vstart,vmin,vmax},DampingFactor->0.8,MaxIterations->100],{i,1,nu}];

vr

N[((vr-RotateRight[vr])/dvj)]

imax=Dimensions[vr][[1]];
vr2=vr^2;

Roots:

imax=;
vr={};
vr2=(vr^2);

Coefficients

sp=1-((a Sin[((vr h)/a)])/(vr h));
cp=1+((a Sin[((ur h)/a)])/(ur h));

Transient Temperature Distribution for Radiant Surface Heating

T[time_,r_,z_]:=(1/(k h))(2 (a^2) (∑_(i=1)^imax)(∑_(m=1)^mmax)(pm[[m]] (1-Exp[-(time k (((ur[[i]])^2)+zetam2[[m]]))/(ro c (a^2))]) Cos[((ur[[i]] h)/(2 a))] Cos[((ur[[i]] z)/a)] BesselJ[0,((zetam[[m]] r)/a)])/((((ur[[i]])^2)+zetam2[[m]]) cp[[i]]))-(1/(k h))(2 (a^2) (∑_(i=1)^imax)(∑_(m=1)^mmax)(pm[[m]] (1-Exp[-(time k (((vr[[i]])^2)+zetam2[[m]]))/(ro c (a^2))]) Sin[((vr[[i]] h)/(2 a))] Sin[((vr[[i]] z)/a)] BesselJ[0,((zetam[[m]] r)/a)])/((((vr[[i]])^2)+zetam2[[m]]) sp[[i]]))

N[T[tc,0,-(h/2)]]

Evolution of the temperature at the center on the surface:

ns=5;
nd=5;
nt=nd ns+1;
to=(1/(10^5));
fact=(10^(1/ns));
N[ts=Table[to (fact^(i-1)),{i,nt}];]
N[temper=Table[T[tc ts[[j]],0,-(h/2)],{j,nt}];]
N[Max[temper]]

ListLogLinearPlot[Transpose[Partition[Join[tc ts,temper],nt]],Joined->True,GridLines->Automatic,PlotStyle->Thickness[0.005]]

Diametrical cuts of the surface coating temperature field at various times:

radDiv=43;
rad=Table[i a/radDiv,{i,0,radDiv}];
N[tempern=Table[T[tc/100,rad[[i]],-(h/2)],{i,1,radDiv+1}];]
N[temper2=Table[T[tc/1000,rad[[i]],-(h/2)],{i,1,radDiv+1}];]
N[temper3=Table[T[tc/10000,rad[[i]],-(h/2)],{i,1,radDiv+1}];]
N[temper4=Table[T[tc/100000,rad[[i]],-(h/2)],{i,1,radDiv+1}];]

tc/100

ListPlot[{Transpose[Partition[Join[rad,tempern],radDiv+1]],
Transpose[Partition[Join[rad,temper2],radDiv+1]],
Transpose[Partition[Join[rad,temper3],radDiv+1]],
Transpose[Partition[Join[rad,temper4],radDiv+1]]},
Joined->True,GridLines->Automatic,PlotRange->All,PlotStyle->Thickness[0.005],PlotLegends->{"tc/100","tc/1,000","tc/10,000","tc/100,000"},Frame->True,FrameLabel->{"radius (m)","HR Surface Temperature (C)"}]

Compares reasonably well with an ANSYS nonlinear, transient thermal reponse for the OCF mirror with 1 Watt absorbed power :

Warning : The balance of this notebook has not been debugged for the OCF mirror reponse.

Wavefront Aberrations due to Index Gradients

Formulation:

c2q[time_,q_,br_]:=(1/(π k h br))((((-1))^q) (2 q+1) 4 P (a^2) dndT (∑_(i=1)^imax)(∑_(m=1)^mmax)((1-Exp[-(time k (ur2[[i]]+zetam2[[m]]))/(ro c (a^2))]) zetam[[m]] Exp[-((zetam2[[m]] (w^2))/(8 (a^2)))] Sin[((ur[[i]] h)/a)] BesselJ[2 q+1,((zetam[[m]] br)/a)])/((zetam2[[m]]+(tau^2)) (zetam2[[m]]+ur2[[i]]) ((BesselJ[0,zetam[[m]]])^2) ur[[i]] cp[[i]]));

psy[time_,r_]:=(1/(k h))(dndT 2 (a^3) (∑_(i=1)^imax)(∑_(m=1)^mmax)(pm[[m]] (1-Exp[-(time k (ur2[[i]]+zetam2[[m]]))/(ro c (a^2))]) Sin[((ur[[i]] h)/a)] BesselJ[0,((zetam[[m]] r)/a)])/((ur2[[i]]+zetam2[[m]]) ur[[i]] cp[[i]]));

Steady-State Aberrations:

Optical Path Distortion:

np=26;
radii=Table[i,{i,0,b,(b/(np-1))}];
N[aber=Table[Abs[psy[10 tc,radii[[i]]]],{i,1,np}];]

aberpp=N[aber[[1]]-aber[[np]]]

N[aber[[1]]]

OPDLensingsurfAbsb=Transpose[Partition[Join[radii,aber],np]];

Optical Path Difference in Thermal Lensing for Surface Absorption (r <= b):

N[OPDLensingsurfAbsb]

np=26;
radii2=Table[i,{i,b,a,((a-b)/(np-1))}];
N[aber2=Table[Abs[psy[10 tc,radii2[[i]]]],{i,1,np}];]
abera=Join[aber,aber2];
radiia=Join[radii,radii2];

aberapp=N[abera[[1]]-abera[[2 np]]]

N[abera[[2 np]]]

OPDLensingsurfAbsa=Transpose[Partition[Join[radiia,abera],2 np]];

Optical Path Difference in Thermal Lensing for Surface Absorption (r <= a):

N[OPDLensingsurfAbsa]

Zernike Expansion:

zernikesb=N[Table[c2q[10 tc,q,b],{q,0,5}]];
MatrixForm[zernikesb]

Table[-(((-1))^i),{i,1,6}].zernikesb

Zb=Table[(∑_(q=0)^5)(zernikesb[[q+1]] Z2q[q,r]/.r->((radii[[i]])/b)),{i,1,np}];

Zfitb=Transpose[Partition[Join[radii,Zb],np]];

zernikesa=N[Table[c2q[10 tc,q,a],{q,0,5}]];
MatrixForm[zernikesa]

Table[-(((-1))^i),{i,1,6}].zernikesa

Za=Table[(∑_(q=0)^5)(zernikesa[[q+1]] Z2q[q,r]/.r->((radiia[[i]])/a)),{i,1,2 np}];

Zfita=Transpose[Partition[Join[radiia,Za],2 np]];

Transient Aberrations:

ns=5;
nd=5;
nt=nd ns+1;
to=(1/(10^5));fact=(10^(1/ns));
N[ts=Table[to (fact^(i-1)),{i,nt}];]

N[Zernike2=Table[Abs[c2q[tc ts[[j]],1,b]],{j,nt}];]
N[Max[Zernike2]]

N[Zernike2=Table[Abs[c2q[tc ts[[j]],1,a]],{j,nt}];]
N[Max[Zernike2]]

N[Zernike10=Table[Abs[c2q[tc ts[[j]],5,b]],{j,nt}];]
N[Max[Zernike10]]

Thermoelastic Analysis

Steady-State Wavefront Aberrations due to Surface Distortion:

The solutions below do not include the St. Venant correction for the effect of the boundary condition in the vicinity of the outer radius since the primary concern is the clear optical region near the optic center, and the St. Venant correction is only about 5%.

Optical Path Distortion:

Sss[r_]:=-(ni (eps (a^2)) Plus@@(pm ((zetam+tau)-(zetam-tau) Exp[-4 gamma]) BesselJ[0,((zetam r)/a)])/(zetam ((((zetam+tau))^2)-(((zetam-tau))^2) Exp[-4 gamma])))/((lambda+mu) 2 k);

np=26;
radii=Table[i,{i,0,b,(b/(np-1))}];
N[surfss=Table[Abs[Sss[radii[[i]]]],{i,1,np}];]

surfsspp=N[surfss[[1]]-surfss[[np]]]

N[surfss[[1]]]

SurfAbsb=Transpose[Partition[Join[radii,surfss-surfss[[np]]],np]];

Surface Displacement for Surface Absorption (r <= b):

N[SurfAbsb]

np=26;
radii2=Table[i,{i,b,a,((a-b)/(np-1))}];
N[surfss2=Table[Abs[Sss[radii2[[i]]]],{i,1,np}];]
surfssa=Join[surfss,surfss2];
radiia=Join[radii,radii2];

surfssapp=N[surfssa[[1]]-surfssa[[2 np]]]

N[surfssa[[2 np]]]

SurfAbsa=Transpose[Partition[Join[radiia,surfssa-surfssa[[2 np]]],2 np]];

Surface Displacement for Surface Absorption (r <= a):

N[SurfAbsa]

Formulation:

tm[time_,z_]:=-(a^2) (∑_(i=1)^imax)((2 (a^2) pm (1-Exp[-(time k (((ur[[i]])^2)+zetam2))/(ro c (a^2))]) Cos[((ur[[i]] h)/(2 a))] Cos[((ur[[i]] z)/a)])/(k h (((((ur[[i]])^2)+zetam2))^2) cp[[i]])-(2 (a^2) pm (1-Exp[-(time k (((vr[[i]])^2)+zetam2))/(ro c (a^2))]) Sin[((vr[[i]] h)/(2 a))] Sin[((vr[[i]] z)/a)])/(k h (((((vr[[i]])^2)+zetam2))^2) sp[[i]]));

tmprime[time_,z_]:=a (∑_(i=1)^imax)((2 ur[[i]] (a^2) pm (1-Exp[-(time k (((ur[[i]])^2)+zetam2))/(ro c (a^2))]) Cos[((ur[[i]] h)/(2 a))] Sin[((ur[[i]] z)/a)])/(k h (((((ur[[i]])^2)+zetam2))^2) cp[[i]])+(2 vr[[i]] (a^2) pm (1-Exp[-(time k (((vr[[i]])^2)+zetam2))/(ro c (a^2))]) Sin[((vr[[i]] h)/(2 a))] Cos[((vr[[i]] z)/a)])/(k h (((((vr[[i]])^2)+zetam2))^2) sp[[i]]));

Zernike Expansion:

N[tmp=tm[10 tc,(h/2)];]
N[tmn=tm[10 tc,-(h/2)];]

N[tmprimep=tmprime[10 tc,(h/2)];]
N[tmprimen=tmprime[10 tc,-(h/2)];]

A1m=-(ni ((tmprimen+tmprimep) Sinh[gamma]+(1/a)(zetam (tmn-tmp) Cosh[gamma])))/((lambda+mu) (gamma-Sinh[gamma] Cosh[gamma]));

A2m=-(ni ((tmprimen-tmprimep) Cosh[gamma]+(1/a)(zetam (tmn+tmp) Sinh[gamma])))/((lambda+mu) (gamma+Sinh[gamma] Cosh[gamma]));

B1m=-(ni ((1/a)(zetam (tmn+tmp) (((mu Sinh[gamma])/(lambda+2 mu))-((lambda+mu) gamma Cosh[gamma])/(lambda+2 mu)))+(tmprimen-tmprimep) (Cosh[gamma]-((lambda+mu) gamma Sinh[gamma])/(lambda+2 mu))))/(2 (lambda+mu) (gamma+Sinh[gamma] Cosh[gamma]));

B2m=-(ni ((1/a)(zetam (tmn-tmp) (((mu Cosh[gamma])/(lambda+2 mu))-((lambda+mu) gamma Sinh[gamma])/(lambda+2 mu)))+(tmprimen+tmprimep) (Sinh[gamma]-((lambda+mu) gamma Cosh[gamma])/(lambda+2 mu))))/(2 (lambda+mu) (gamma-Sinh[gamma] Cosh[gamma]));

Bm=N[-B1m Sinh[((zetam (-h))/(2 a))]-B2m Cosh[((zetam (-h))/(2 a))]+((ni tmprimen)/(lambda+2 mu))+((lambda+3 mu) (A1m Cosh[((zetam (-h))/(2 a))]+A2m Sinh[((zetam (-h))/(2 a))]))/(2 (lambda+2 mu))-((lambda+mu) zetam (-h) (A1m Sinh[((zetam (-h))/(2 a))]+A2m Cosh[((zetam (-h))/(2 a))]))/(2 (lambda+2 mu) 2 a)];

ct2q[q_,br_]:=(((-1))^q) (4 q+2) a (∑_(m=1)^mmax)(Bm[[m]] BesselJ[2 q+1,((zetam[[m]] br)/a)])/(zetam[[m]] br)

Cssb=N[Table[Abs[ct2q[q,b]],{q,0,5}]];
MatrixForm[Cssb]

Table[-(((-1))^i),{i,1,6}].Cssb

Zb=Table[(∑_(q=0)^5)(Cssb[[q+1]] Z2q[q,r]/.r->((radii[[i]])/b)),{i,1,np}];

Zfitb=Transpose[Partition[Join[radii,Zb],np]];

Cssa=N[Table[ct2q[q,a],{q,0,5}]];
MatrixForm[Cssa]

Table[-(((-1))^i),{i,1,6}].Cssa

Za=Table[(∑_(q=0)^5)(Cssa[[q+1]] Z2q[q,r]/.r->((radiia[[i]])/a)),{i,1,2 np}];

Zfita=Transpose[Partition[Join[radiia,-(Za/2)],2 np]];

Transient Wavefront Aberrations due to Surface Distortion

TBD

Optical Path Distortion:

xip=((ur h)/(2 a));
etap=((vr h)/(2 a));

Spm[i_,m_,time_]:=(2 (a^2) pm[[m]] (1-Exp[-(time k (((ur[[i]])^2)+zetam2[[m]]))/(ro c (a^2))]) Cos[((ur[[i]] h)/(2 a))])/(k h (((ur[[i]])^2)+zetam2[[m]]) cp[[i]]);
Apm[i_,m_,time_]:=-(2 (a^2) pm[[m]] (1-Exp[-(time k (((vr[[i]])^2)+zetam2[[m]]))/(ro c (a^2))]) Sin[((vr[[i]] h)/(2 a))])/(k h (((vr[[i]])^2)+zetam2[[m]]) sp[[i]]);

S[time_,r_]:=-(a ni (∑_(i=1)^imax)(∑_(m=1)^mmax)((Sinh[gamma[[m]]] Spm[i,m,time] (ur[[i]] Cosh[gamma[[m]]] Sin[etap[[i]]]+zetam[[m]] Sinh[gamma[[m]]] Cos[xip[[i]]]))/((Sinh[gamma[[m]]] Cosh[gamma[[m]]]+gamma[[m]]) (((ur[[i]])^2)+zetam2[[m]]))+(Cosh[gamma[[m]]] Apm[i,m,time] (vr[[i]] Sinh[gamma[[m]]] Cos[etap[[i]]]-zetam[[m]] Cosh[gamma[[m]]] Sin[etap[[i]]]))/((Sinh[gamma[[m]]] Cosh[gamma[[m]]]-gamma[[m]]) (((vr[[i]])^2)+zetam2[[m]]))) BesselJ[0,((zetam[[m]] r)/a)])/(lambda+mu);

N[S[10 tc,0]]

N[S[10 tc,a]]

np=26;
radii=Table[i,{i,0,b,(b/(np-1))}];
N[Sss=Table[Abs[S[10 tc,radii[[i]]]],{i,1,np}];]

Sss_pp=N[Sss[[1]]-Sss[[np]]]

Bulk Absorption

Steady-State Thermal Field Solution

T[r_,z_]:=(1/k)(alpha (a^2) (∑_(m=1)^mmax)(pm[[m]] (1-(2 tau Cosh[((zetam[[m]] z)/a)])/((zetam[[m]]+tau) Exp[((zetam[[m]] h)/(2 a))]-(zetam[[m]]-tau) Exp[-((zetam[[m]] h)/(2 a))])) BesselJ[0,((zetam[[m]] r)/a)])/zetam2[[m]])

N[T[0,-(h/2)]]

temp=Table[T[i,j],{i,-a,a,(a/24)},{j,-(h/2),(h/2),(h/16)}];

Max[temp]

Mean[Flatten[temp]]

Transient Thermal Field Solution

T[time_,r_,z_]:=(∑_(i=1)^imax)(∑_(m=1)^mmax)(4 alpha (a^3) pm[[m]] (((Sin[((ur[[i]] h)/(2 a))])/(ur[[i]]))-((tau Cos[((ur[[i]] h)/(2 a))])/(zetam2[[m]]+ur2[[i]]))) (1-Exp[-(k time (zetam2[[m]]+ur2[[i]]))/(ro c (a^2))]) Cos[((ur[[i]] z)/a)] BesselJ[0,((zetam[[m]] r)/a)])/(k h zetam2[[m]] cp[[i]])

N[T[10 tc,0,0]]

ns=5;
nd=5;
nt=nd ns+1;
to=(1/(10^5));
fact=(10^(1/ns));
ts=N[Table[to (fact^(i-1)),{i,nt}]];

N[temper=Table[T[tc ts[[j]],0,-(h/2)],{j,nt}];]
N[Max[temper]]

TransientTemp=Transpose[Partition[Join[(tc ts)/3600,temper],nt]];

Evolution of the temperature at the center on the surface:

Diametrical cuts of the coating temperature field at various times:

Axial cuts of the temperature field at various times:

N[temper4=Table[T[tc,0,j],{j,-(h/2),(h/2),(h/48)}];]
N[temper41=Table[T[(tc/1000),0,j],{j,-(h/2),(h/2),(h/48)}];]

Thermal Aberrations due to Index Gradients

Formulation:

c2q[time_,q_,br_]:=(1/(π k (h^2) br))((((-1))^q) (2 q+1) 8 P (a^3) dndT (∑_(i=1)^imax)(∑_(m=1)^mmax)((1-Exp[-(time k (ur2[[i]]+zetam2[[m]]))/(ro c (a^2))]) Exp[-((zetam2[[m]] (w^2))/(8 (a^2)))] (((1-Cos[((ur[[i]] h)/a)])/(ur[[i]]))-((tau Sin[((ur[[i]] h)/a)])/(ur2[[i]]+zetam2[[m]]))) BesselJ[2 q+1,((zetam[[m]] br)/a)])/((zetam2[[m]]+(tau^2)) zetam[[m]] ur[[i]] ((BesselJ[0,zetam[[m]]])^2) cp[[i]]));

psy[time_,r_]:=(1/(k h))(dndT 4 alpha (a^4) (∑_(i=1)^imax)(∑_(m=1)^mmax)(pm[[m]] (1-Exp[-(time k (ur2[[i]]+zetam2[[m]]))/(ro c (a^2))]) (((1-Cos[((ur[[i]] h)/a)])/(ur[[i]]))-((tau Sin[((ur[[i]] h)/a)])/(ur2[[i]]+zetam2[[m]]))) BesselJ[0,((zetam[[m]] r)/a)])/(zetam2[[m]] cp[[i]] ur[[i]]))

Steady-State Aberrations:

Optical Path Distortion:

np=26;
radii=Table[i,{i,0,b,(b/(np-1))}];
N[bulkaber=Table[Abs[psy[10 tc,radii[[i]]]],{i,1,np}];]
bulkaberpp=N[bulkaber[[1]]-bulkaber[[np]]]

N[bulkaber[[1]]]

N[bulkaber[[np]]]

OPDBulkAbsb=Transpose[Partition[Join[radii,bulkaber],np]];

np=26;
radii2=Table[i,{i,b,a,((a-b)/(np-1))}];
N[bulkaber2=Table[Abs[psy[10 tc,radii2[[i]]]],{i,1,np}];]
bulkabera=Join[bulkaber,bulkaber2];
radiia=Join[radii,radii2];
bulkaberapp=N[bulkabera[[1]]-bulkabera[[2 np]]]

N[bulkabera[[2 np]]]

OPDBulkAbsa=Transpose[Partition[Join[radiia,bulkabera],2 np]];

Reflection OPD for Bulk Absorption (r <= a):

N[OPDBulkAbsa]

Zernike Expansion:

zernikesb=N[Table[c2q[10 tc,q,b],{q,0,5}]];
MatrixForm[zernikesb]

Table[-(((-1))^i),{i,1,6}].zernikesb

zernikesa=N[Table[c2q[10 tc,q,a],{q,0,5}]];
MatrixForm[zernikesa]

Table[-(((-1))^i),{i,1,6}].zernikesa

Transient Aberrations:

ns=5;
nd=5;
nt=nd ns+1;
to=(1/(10^5));
fact=(10^(1/ns));
N[ts=Table[to (fact^(i-1)),{i,nt}];]

N[Zernike2=Table[Abs[c2q[tc ts[[j]],1]],{j,nt}];]
N[Max[Zernike2]]

N[Zernike10=Table[Abs[c2q[tc ts[[j]],5]],{j,nt}];]
N[Max[Zernike10]]

Thermoelastic Analysis

Steady State Aberrations

Optical Path Distortion:

Sss[r_]:=-(a ni (∑_(i=1)^imax)(∑_(m=1)^mmax)(Sinh[gamma[[m]]] 4 (a^3) alpha pm[[m]] (((Sin[((ur[[i]] h)/(2 a))])/(ur[[i]]))-((tau Cos[((ur[[i]] h)/(2 a))])/(ur2[[i]]+zetam2[[m]]))) (ur[[i]] Cosh[gamma[[m]]] Sin[((ur[[i]] h)/(2 a))]+zetam[[m]] Sinh[gamma[[m]]] Cos[((ur[[i]] h)/(2 a))]) BesselJ[0,((zetam[[m]] r)/a)])/((Sinh[gamma[[m]]] Cosh[gamma[[m]]]+gamma[[m]]) k h (cp[[i]] zetam2[[m]] (ur2[[i]]+zetam2[[m]]))))/(lambda+mu)

np=26;
radii=Table[i,{i,0,b,(b/(np-1))}];
N[bulkaber=Table[Abs[Sss[radii[[i]]]],{i,1,np}];]
bulkaberpp=N[bulkaber[[1]]-bulkaber[[np]]]

N[bulkaber[[1]]]

N[bulkaber[[np]]]

SurfBulkAbsb=Transpose[Partition[Join[radii,bulkaber],np]];

Surface Displacement for Bulk Absorption (r <= b):

N[SurfBulkAbsb]

np=26;
radii2=Table[i,{i,b,a,((a-b)/(np-1))}];
N[bulkaber2=Table[Abs[Sss[radii2[[i]]]],{i,1,np}];]
bulkssa=Join[bulkaber,bulkaber2];
radiia=Join[radii,radii2];

bulkssapp=N[bulkssa[[1]]-bulkssa[[2 np]]]

N[bulkssa[[2 np]]]

SurfBulkAbsa=Transpose[Partition[Join[radiia,bulkssa-bulkssa[[2 np]]],2 np]];

Surface Displacement for Bulk Absorption (r <= a):

N[SurfBulkAbsa]

Formulation:

tm[time_,z_]:=-(a^2) (∑_(i=1)^imax)(4 alpha (a^3) pm (1-Exp[-((time k (ur2[[i]]+zetam2))/(ro c (a^2)))]) (-((tau Cos[((ur[[i]] h)/(2 a))])/(ur2[[i]]+zetam2))+((Sin[((ur[[i]] h)/(2 a))])/(ur[[i]]))) Cos[((ur[[i]] z)/a)])/(k h zetam2 cp[[i]] (ur2[[i]]+zetam2));

tmprime[time_,z_]:=a (∑_(i=1)^imax)(4 alpha (a^3) pm (1-Exp[-((time k (ur2[[i]]+zetam2))/(ro c (a^2)))]) (-((tau Cos[((ur[[i]] h)/(2 a))])/(ur2[[i]]+zetam2))+((Sin[((ur[[i]] h)/(2 a))])/(ur[[i]]))) Sin[((ur[[i]] z)/a)])/(k h zetam2 cp[[i]] (ur2[[i]]+zetam2));

Zernike Expansion:

N[tmp=tm[10 tc,(h/2)];]
N[tmn=tm[10 tc,-(h/2)];]

N[tmprimep=tmprime[10 tc,(h/2)];]
N[tmprimen=tmprime[10 tc,-(h/2)];]

A1m=-(ni ((tmprimen+tmprimep) Sinh[gamma]+(1/a)(zetam (tmn-tmp) Cosh[gamma])))/((lambda+mu) (gamma-Sinh[gamma] Cosh[gamma]));

A2m=-(ni ((tmprimen-tmprimep) Cosh[gamma]+(1/a)(zetam (tmn+tmp) Sinh[gamma])))/((lambda+mu) (gamma+Sinh[gamma] Cosh[gamma]));

B1m=-(ni ((1/a)(zetam (tmn+tmp) (((mu Sinh[gamma])/(lambda+2 mu))-((lambda+mu) gamma Cosh[gamma])/(lambda+2 mu)))+(tmprimen-tmprimep) (Cosh[gamma]-((lambda+mu) gamma Sinh[gamma])/(lambda+2 mu))))/(2 (lambda+mu) (gamma+Sinh[gamma] Cosh[gamma]));

B2m=-(ni ((1/a)(zetam (tmn-tmp) (((mu Cosh[gamma])/(lambda+2 mu))-((lambda+mu) gamma Sinh[gamma])/(lambda+2 mu)))+(tmprimen+tmprimep) (Sinh[gamma]-((lambda+mu) gamma Cosh[gamma])/(lambda+2 mu))))/(2 (lambda+mu) (gamma-Sinh[gamma] Cosh[gamma]));

Bm=N[-B1m Sinh[((zetam (-h))/(2 a))]-B2m Cosh[((zetam (-h))/(2 a))]+((ni tmprimen)/(lambda+2 mu))+((lambda+3 mu) (A1m Cosh[((zetam (-h))/(2 a))]+A2m Sinh[((zetam (-h))/(2 a))]))/(2 (lambda+2 mu))-((lambda+mu) zetam (-h) (A1m Sinh[((zetam (-h))/(2 a))]+A2m Cosh[((zetam (-h))/(2 a))]))/(2 (lambda+2 mu) 2 a)];

ct2q[q_,br_]:=(((-1))^q) (4 q+2) a (∑_(m=1)^mmax)(Bm[[m]] BesselJ[2 q+1,((zetam[[m]] br)/a)])/(zetam[[m]] br)

Cssb=N[Table[ct2q[q,b],{q,0,5}]];
MatrixForm[Cssb]

Table[-(((-1))^i),{i,1,6}].Cssb

Cssa=N[Table[ct2q[q,a],{q,0,5}]];
MatrixForm[Cssa]

Table[-(((-1))^i),{i,1,6}].Cssa

General Transient Solution

S[time_,r_]:=-(a ni (∑_(i=1)^imax)(∑_(m=1)^mmax)(Sinh[gamma[[m]]] 4 (a^3) alpha pm[[m]] (1-Exp[-(time k (ur2[[i]]+zetam2[[m]]))/(ro c (a^2))]) (((Sin[((ur[[i]] h)/(2 a))])/(ur[[i]]))-((tau Cos[((ur[[i]] h)/(2 a))])/(ur2[[i]]+zetam2[[m]]))) (ur[[i]] Cosh[gamma[[m]]] Sin[((ur[[i]] h)/(2 a))]+zetam[[m]] Sinh[gamma[[m]]] Cos[((ur[[i]] h)/(2 a))]) BesselJ[0,((zetam[[m]] r)/a)])/((Sinh[gamma[[m]]] Cosh[gamma[[m]]]+gamma[[m]]) k h (cp[[i]] zetam2[[m]] (ur2[[i]]+zetam2[[m]]))))/(lambda+mu)

N[S[10 tc,0]]

Comparisons

Surface Absorption:
Thermal Lensing (Transmission) vs. Thermoelastic Distortion (Reflection)

OPDSurfAbsa=Transpose[Partition[Join[radiia,2 (surfssa-surfssa[[2 np]])],2 np]];

Bulk Absorption:
Thermal Lensing (Transmission) vs. Thermoelastic Distortion (Reflection)

OPDReflBulkAbsa=Transpose[Partition[Join[radiia,2 (bulkssa-bulkssa[[2 np]])],2 np]];

Thermal Lensing (Transmission):
Surface vs. Bulk Absorption

Thermoelastic Distortion (Reflection):
Surface vs. Bulk Absorption