% TE calculation using 2D axisymmetric model. This script calls the model % file, extracts grad_T and does the calculations for the different % frequency % All parameter are defined in this script and passed to the COMSOL model % via arguments % Parameters, Quantitites in [m] param.beam.radius = 0.05; param.mirror.radius = 100 * param.beam.radius; %Define the radius of the mirror as 100 times that of the beam radius param.mirror.height = 200 * param.beam.radius; % Define the height of the cylinder to be 200 times the beam radius param.COMSOL.total_radial_slices = 20; param.COMSOL.total_z_slices = 40; param.constant.boltzmann = 1.38065e-23; param.constant.F0 = 1e-2;% Force amplitude param.constant.T_amb = 300; param.COMSOL.freq = 10.^[0:0.25:4];% frequency range param.COMSOL.simulation_cycles = 10;% The number of cycles of Pressure Osc. param.COMSOL.simulation_time_slice_per_cycle = 20; % Time steps in each cycle % Parameters for fused silica param.material = fused_silicadata; % Assign values to different varibales alpha = param.material.alpha; E = param.material.E; C = param.material.specific_heat; kappa = param.material.thermal_conductivity; sigma = param.material.sigma; rho = param.material.density; kB = param.constant.boltzmann; T = param.constant.T_amb; r0 = param.beam.radius; F0 = param.constant.F0; disp('Starting Time'); tic; freq = param.COMSOL.freq; psd = zeros(1,length(freq));% to store the PSD values %%Calculation and analysis : for i = 1:length(freq) param.COMSOL.fmod = freq(i);% COMSOL will be asked to simulate with this frequency dt = 1/(freq(i) * param.COMSOL.simulation_time_slice_per_cycle); param.COMSOL.simulation_time_step = dt;% time steps param.COMSOL.simulation_end_time = 1/freq(i) * param.COMSOL.simulation_cycles;% run simulation until this time t_end = param.COMSOL.simulation_end_time; disp(['Running COMSOL for frequency = ' num2str(freq(i))]); %% Call to COMSOL model = almost_infinite_2D_axisym_TE_4(param); %% Data Extraction % definitions r_slices = param.COMSOL.total_radial_slices;% assign new variables for convinience dR = param.mirror.radius/r_slices;% size of each R slice z_slices = param.COMSOL.total_z_slices; dz = param.mirror.height/z_slices;% size of each z slice t_slices = t_end/dt; disp('## Extracting data'); % Store the data for each r, z for the last three cycles of simulation % time dat = zeros(r_slices,z_slices,floor(3/(freq(i)*dt)));% to store data from COMSOL for each r, z and 't'- for the last three cycles for m = 0:(r_slices-1) for k = 0:(z_slices-1) dat(m+1,k+1,:) = mphinterp(model,'ht.gradTmag','coord',[m*dR ; k*dz],'t', [(t_slices - floor(3/(freq(i)*dt)) + 1) : t_slices]); % The data is extracted for the last 3 cycles of simulation end end t_eval = dt * ((t_slices - floor(3/(freq(i)*dt)) + 1) : t_slices);% the time values for the last 3 cycles; temp = exp(j * 2*pi*freq(i) .* t_eval);% stores all exp(i * omega * t) values kernel = zeros(1,1,length(temp)); kernel(1,1,:) = temp;% Now dimensions of kernel and dat_fc matches % Without the above modification, the assignment in the nested for % gives an error dat_fc = zeros(r_slices,z_slices);% matrix to store the fourier coefficients for each r, z time signal %% Fourier coefficients corresponding to fmod are extracted disp('## Extracting Fourier coefficients'); for m = 1:r_slices for k = 1:z_slices dat_fc(m,k) = abs(2*mean(kernel.*dat(m,k,:)));% this is sqrt(an^2 + bn^2) end end disp('## Extraction Complete. Evaluating volume integral'); %% Calculation of the integral of grad_T^2 dat_fc = dat_fc.^2; % square the values since we need (an^2 + bn^2) instead of its sqrt dat_fc = dat_fc';% take a transpose so that going down a column you get all z values for a fixed r z = 0:dz:(param.mirror.height-dz);% the z values r = 0:dR:(param.mirror.radius-dR);% the r values z_integral = trapz(z,dat_fc);% integrate over the columns of the matrix dat_fc vol_integral = trapz(r,2*pi*r.*z_integral);% total volume integration of (an^2 + bn^2) % There must be an additional factor of 1/2 since the integral we want % is 1/2 (an^2 + bn^2) vol_integral = 0.5 * vol_integral;% Now we have time_avg as calculated in the paper psd(i) = vol_integral/(2*pi*freq(i))^2;% divide by omega^2 and store toc; end psd = 8*kB*kappa/F0^2 * psd; %Multiply all the prefactors in the formula of FDT loglog(freq, psd,'b.'); hold on; LT_constant = (8*(1 + sigma)^2*kappa*alpha^2*kB*T^2) / (sqrt(2*pi)*(rho*C)^2*r0^3); loglog(freq, LT_constant./(2*pi.*freq).^2,'ro'); hold off;