ID 
Date 
Author 
Type 
Category 
Subject 
39

Thu Jun 27 17:12:33 2013 
Arnaldo Rodriguez  Optics  General  Defocusing in PControlled Manual Simulation 
To test the effectiveness of the PID loop, I've only given a nonzero value to the K_p constant and varied it over a series of runs to verify the response by the controller.
The following is a plot of the baseline compared to the Pactive simulations for different proportionality gain constants. I've yet to determine whether all the trends in behavior are physical.
Each different case has an evident droop (steadystate error), and seem to be destabilizing as K_p increases. However, most cases reach steadystate quicker than the baseline.
It should be noted that the ring heater is not correcting any effects from the beam, as it is turned off. The ring heater is correcting itself; or rather, its own initial heat pulse of 5 W * delta_t (600s) = 3000 J before the controller becomes active. 
40

Fri Jun 28 11:54:10 2013 
Emory Brown  Optics  General  Maintaining a constant volume while modifying the radii of the test mass faces 
When changing the dimensions of the mirror, we want to keep a constant volume. To do so, we must modify the height of the test mass. The area of the test mass if R is the radius of the front mirror face and r the radius of the back mirror face is (1/3)*pi*h*(R^2+R*r+r^2)=area, so using the values currently in our model, R=r=0.17 m and h=0.2 m, area=0.0181584055 m^3=pi*0.2*0.17^2 m^3, so to ensure that the volume of the test mass remains constant, h=0.01734/(R^2+R*r+r^2). Note that this does give 0.2 m when R=r=0.17 m. 
41

Fri Jun 28 14:03:07 2013 
Deep Chatterjee  General  General  Approach for TE noise for beam in cavity 
In this post I mention of the approach that would be taken to solve for the TE noise spectrum for a beam inside the substrate cavity.
Completing the same would be my plan for the next week.
In the COMSOL simulation, a gaussian pressure will be applied on both the flat surface of the cylinder. The strain would be calculated
using the Structural mechanics module and the temperature field would be calculated using the Heat Transfer module. The approach
is similar to that followed in Liu and Thorne. The Work dissipated would be calculated using the gradient of T. The spectral density is
then calculated using the FDT.
In the attachment, I have written briefly about the material I have studied this week along with my thoughts on the problem of the TE noise
calculation. 
Attachment 1: LL_derivations.pdf


42

Mon Jul 1 13:01:58 2013 
Arnaldo Rodriguez  Optics  General  Brief Parameter Study on Pactive Control Loop 
The following is a plot of the defocus as a function of time with only the proportional gain being nonzero. Only three cases were run before the program crashed.
The dotted line indicates the setpoint.
One can clearly see the loop becomes unstable and oscillates out of control at a proportional gain of 240, and can see the characteristic steadystate offset of the Pcontroller from the set point develop in the other cases. 
43

Mon Jul 1 17:31:40 2013 
Arnaldo Rodriguez  Optics  General  Parameter Study on Sampling Time for Pactive Loop 
The following is a plot of the defocus as a function of time for a fixed proportionality gain of 150, and different listed sampling times.
The discrepancy between the behaviors is enough to warrant a proper mathematical explanation; the max. overshoot, for example, is inversely proportional to the sampling rate.

44

Tue Jul 2 13:32:30 2013 
Emory Brown  Optics  General  Using Matlab to run our COMSOL simulations 
We want to be able to test a number of values for the ratio of the radii of the frustum which will be much easier once we can get Matlab to automate running our COMSOL model. I have been able to output the COMSOL model to a .m file and was able to figure out a way to export the Umax value we want to Matlab (see the attached .m file).
The model ran in Matlab after debugging it for a short while and returned the same result as the COMSOL model had by itself. I setup a script to loop over some values of RatioRadii and feed them to the model. The first value it fed to the model worked and returned 1.2099e10 J for RatioRadii=0.1, but for RatioRadii=0.2 the model had a
meshing error. I am going to have to see if there is some way we can modify the mesh to fix this error. When I briefly looked into it today, I found that for several RatioRadii which I tried, the object could be meshed for certain finenesses, but not others. I am going to spend tomorrow working on the first progress report, but after that is done I will
work on fixing this meshing error. 
Attachment 1: model_7_1_13.m

function [Umax] = model(ratioRadii)
% ratioRadii must be a scalar
% model_7_1_13.m
%
% Model exported on Jul 1 2013, 17:08 by COMSOL 4.3.2.189.
import com.comsol.model.*
import com.comsol.model.util.*
model = ModelUtil.create('Model');
... 220 more lines ...

Attachment 2: looper.m

function [RatioRValues,UmaxValues] = looper
% Iterates over a range of RatioRValues
RatioRValues=[0.1:0.1:10]
UmaxValues=[]
for r=RatioRValues,
UmaxValues = [UmaxValues model_7_1_13(r)];
UmaxValues
end
... 4 more lines ...

45

Wed Jul 3 00:04:46 2013 
Koji  Optics  General  Using Matlab to run our COMSOL simulations 
Please use zipped files instead of uploading extremely long text files on the elog.
Or use wiki or svn instead of elog for uploading them.
And appropriate word wrapping so that we can read the entries. Please. 
46

Thu Jul 4 15:23:47 2013 
Deep Chatterjee  Optics  General  Attempt to compute TE noise following Liu and Thorne 
Liu and Thorne, in their paper have analytically solved for the TE noise spectra considering the test mass to be a half infinite space.
An attempt was made to get similar results using a COMSOL model. In this post, the rough profile of the PSD plot is given and is compared
to Liu and Thorne's result.
Steps followed by Liu and Thorne:
> Oscillating pressure, scaled according to beam profile, is applied one of the surfaces.
> The stress balance equation is solved borrowing algorithm from Landau Lifshitz.
> The stress and the temperature perturbation is related using an equation given in Chapter 1 of LL
> The W_diss is calculated from the gradient of the temperature.
> PSD is calculated using FDT due to Levin
> The resulting plot looks like the following  ***Only the frequency dependence is plotted here leaving out the prefactors which will
appear just as an additive constant in the figure. The aim was to see how the plot of the simulation looks like.
Steps followed in the COMSOL model:
> 2D axisymmetric model was used since the applied was applied radially symmetrically.
> The dimensions of the test mass(cylinder) was defined larger as compared to the beam spot radius.
> The physics added was Solid Stress and Heat Transfer in Solids.
> Thermal expansion was enabled in the Solid Stress physics
> Pressure work was enabled in Heat Transfer physics which takes into consideration heat due to stresses i.e. adds this heat source to the heat eqn
> Heat insulation wsa applied to all surfaces.
> All surfaces apart from the one where pressure was applied was held fixed to mimic infinity
> A custom mesh was put.
> The model was exported as '.m'
> It was called for frequency range 10^(0 : 0.25 : 4)
> Following Liu and Thorne, the temperature gradient was integrated over the volume using mphint2() function in Matlab
> The time average was taken over the second half of the cycle to avoid transient which are assumed to die off after a couple of initial cycles.
> All the prefactors were ignored while plotting
> The plot shown below is integral{grad_T}^2} / omega^2 which is proportional to PSD
> We just wanted to see if at all we got anything looking similar to the analytic result of Liu and Thorne (shown above)
The above plot is when the radius is 100 times the beam spot and height of the cylinder is 200 times the beam spot
Other values of radius and height were also tried. The one below shows similar plot as above but for radius 10 times
the beam spot and height 20 times the same
While the one below is for radius = 1000 * beam spot and height = 2000 * beam spot
Here, it is seen that the straight line with a negative slope is seen to some extent.
However, the prefactors, when put in, will show the actual PSD.
Any comments or suggestions related to the model construction in COMSOL or otherwise is welcome. 
Attachment 1: Liu&Thorne_type.bmp

Attachment 5: 2D_almost_infinte_TE.mph

47

Fri Jul 5 14:06:24 2013 
Deep Chatterjee  Optics  General  Attempt to compute TE noise following Liu and Thorne 
Quote: 
Liu and Thorne, in their paper have analytically solved for the TE noise spectra considering the test mass to be a half infinite space.
An attempt was made to get similar results using a COMSOL model. In this post, the rough profile of the PSD plot is given and is compared
to Liu and Thorne's result.
Steps followed by Liu and Thorne:
> Oscillating pressure, scaled according to beam profile, is applied one of the surfaces.
> The stress balance equation is solved borrowing algorithm from Landau Lifshitz.
> The stress and the temperature perturbation is related using an equation given in Chapter 1 of LL
> The W_diss is calculated from the gradient of the temperature.
> PSD is calculated using FDT due to Levin
> The resulting plot looks like the following  ***Only the frequency dependence is plotted here leaving out the prefactors which will
appear just as an additive constant in the figure. The aim was to see how the plot of the simulation looks like.
Steps followed in the COMSOL model:
> 2D axisymmetric model was used since the applied was applied radially symmetrically.
> The dimensions of the test mass(cylinder) was defined larger as compared to the beam spot radius.
> The physics added was Solid Stress and Heat Transfer in Solids.
> Thermal expansion was enabled in the Solid Stress physics
> Pressure work was enabled in Heat Transfer physics which takes into consideration heat due to stresses i.e. adds this heat source to the heat eqn
> Heat insulation wsa applied to all surfaces.
> All surfaces apart from the one where pressure was applied was held fixed to mimic infinity
> A custom mesh was put.
> The model was exported as '.m'
> It was called for frequency range 10^(0 : 0.25 : 4)
> Following Liu and Thorne, the temperature gradient was integrated over the volume using mphint2() function in Matlab
> The time average was taken over the second half of the cycle to avoid transient which are assumed to die off after a couple of initial cycles.
> All the prefactors were ignored while plotting
> The plot shown below is integral{grad_T}^2} / omega^2 which is proportional to PSD
> We just wanted to see if at all we got anything looking similar to the analytic result of Liu and Thorne (shown above)
The above plot is when the radius is 100 times the beam spot and height of the cylinder is 200 times the beam spot
Other values of radius and height were also tried. The one below shows similar plot as above but for radius 10 times
the beam spot and height 20 times the same
While the one below is for radius = 1000 * beam spot and height = 2000 * beam spot
Here, it is seen that the straight line with a negative slope is seen to some extent.
However, the prefactors, when put in, will show the actual PSD.
Any comments or suggestions related to the model construction in COMSOL or otherwise is welcome.

Here is a plot showing the profiles of Liu and Thorne and the COMSOL simulation. Since the constants were not taken into account, one expects the plot to be shifted by a constant equal to the logarithm of the prefactor that has been ignored.
The green 'o' is just a 1/omega^2 plot in the loglog scale which is how Liu and Thorne's result go (this is not linearized PSD, but is doesn't matter as long as we want to see the profiles). The other plot in blue is the COMSOL simulation.
The plot below is a difference between them. One expects it to be a constant but there are a few jitters. Once again, it is emphasized that the numerical values do not give any information related to the actual difference which will only be calculable
once the correct physical constants are put in the formula. I will post another entry once that is done

Attachment 2: L&T_COMSOL_difference.bmp

48

Mon Jul 8 15:47:23 2013 
Emory Brown  Optics  General  Meshing error when varying geometry 
I attempted to run the Matlab code built last week, but received an error that the mesh failed to generate on the inner domain/failed to respect edge element on geometry edge. The error occurred on the second computation with ratioR=0.2, but for ratioR=0.1 the simulation completed successfully. I reproduced the error in COMSOL and noted that for meshes with size at least as course as "course" or extremely fine, a mesh could be constructed, but for sizes in between an error was returned. This region of sizes which throw errors varies for different values of ratioR. By specifying an edge selection of the central line running down the frustum and another along the top edge of the inner cylinder and meshing points on those edges first the rest of the meshing is forced to use those nodes as a seed.
I also noted our COMSOL model had an error: ratioR was applied to the wrong face. This was corrected as was the .m file.
I added an eigenfrequency study and ran it for several values of ratioR, obtaining the following:
ratioR Lowest Drumhead mode (Hz) Lowest Butterfly mode (Hz)
0.1 3624.8 3802.4
0.3 4375.9 3785.4
0.5 4366.5 3650.0
.75 3709.0 3041.85
1 3000 2186.3
1.1 2736.2 1893.0
1.3 2247.6 1426.4
The mirror's front face is 0.17 [m]. Its Gaussian beam size is 0.0156 [m]. Its height is given by 0.01734[m^3]/(R^2*(1+ratioR+ratioR^2)) such that when the two faces have the same radii it is 0.2 [m] long. 
49

Mon Jul 8 18:39:31 2013 
Matt A.  Optics  General  Meshing error when varying geometry 
Quote: 
I attempted to run the Matlab code built last week, but received an error that the mesh failed to generate on the inner domain/failed to respect edge element on geometry edge. The error occurred on the second computation with ratioR=0.2, but for ratioR=0.1 the simulation completed successfully. I reproduced the error in COMSOL and noted that for meshes with size at least as course as "course" or extremely fine, a mesh could be constructed, but for sizes in between an error was returned. This region of sizes which throw errors varies for different values of ratioR. By specifying an edge selection of the central line running down the frustum and another along the top edge of the inner cylinder and meshing points on those edges first the rest of the meshing is forced to use those nodes as a seed.
I also noted our COMSOL model had an error: ratioR was applied to the wrong face. This was corrected as was the .m file.
I added an eigenfrequency study and ran it for several values of ratioR, obtaining the following:
ratioR Lowest Drumhead mode (Hz) Lowest Butterfly mode (Hz)
0.1 3624.8 3802.4
0.3 4375.9 3785.4
0.5 4366.5 3650.0
.75 3709.0 3041.85
1 3000 2186.3
1.1 2736.2 1893.0
1.3 2247.6 1426.4
The mirror's front face is 0.17 [m]. Its Gaussian beam size is 0.0156 [m]. Its height is given by 0.01734[m^3]/(R^2*(1+ratioR+ratioR^2)) such that when the two faces have the same radii it is 0.2 [m] long.

Great Emory,
Can you put these in a plot so it's easier to read? 
50

Mon Jul 8 19:04:35 2013 
Matt A.  Optics  General  Meshing error when varying geometry 
Quote: 
I attempted to run the Matlab code built last week, but received an error that the mesh failed to generate on the inner domain/failed to respect edge element on geometry edge. The error occurred on the second computation with ratioR=0.2, but for ratioR=0.1 the simulation completed successfully. I reproduced the error in COMSOL and noted that for meshes with size at least as course as "course" or extremely fine, a mesh could be constructed, but for sizes in between an error was returned. This region of sizes which throw errors varies for different values of ratioR. By specifying an edge selection of the central line running down the frustum and another along the top edge of the inner cylinder and meshing points on those edges first the rest of the meshing is forced to use those nodes as a seed.
I also noted our COMSOL model had an error: ratioR was applied to the wrong face. This was corrected as was the .m file.
I added an eigenfrequency study and ran it for several values of ratioR, obtaining the following:
ratioR Lowest Drumhead mode (Hz) Lowest Butterfly mode (Hz)
0.1 3624.8 3802.4
0.3 4375.9 3785.4
0.5 4366.5 3650.0
.75 3709.0 3041.85
1 3000 2186.3
1.1 2736.2 1893.0
1.3 2247.6 1426.4
The mirror's front face is 0.17 [m]. Its Gaussian beam size is 0.0156 [m]. Its height is given by 0.01734[m^3]/(R^2*(1+ratioR+ratioR^2)) such that when the two faces have the same radii it is 0.2 [m] long.

When we talked today, I thought that these frequencies were a little low, so I ran a really simple COMSOL model of a cylinder with radius 0.17 m and height 0.2 m (attached).
I solved for the first 10 modes (the first 6 being trivial), and the frequencies I got are listed below:
Mode # Freq [Hz] Type
7 5950 Butterfly
8 5950 Butterfly
9 8106 Drumhead
10 8226 Shear
These are more in the frequency range that I expected.
**What material are you using? I used fused silica.
**I used 'finer' meshing. Perhaps your more complicated meshing is screwing with the solutions, although I'd expect it to make it stiffer, not softer.
**Are you applying force during the eigenvalue calculation? That might make it softer.
**Perhaps some of your constraints are moving the modes around. My model has zero constraints: it's just a cylinder floating in space.
Lets talk about this tomorrow.

Attachment 1: SimpleCylinder.mph

51

Tue Jul 9 14:45:04 2013 
Emory Brown  Optics  General  Meshing error when varying geometry 
I tried removing the constraints and applied force and was able to get close to the same frequencies for the butterfly modes as the simple model. By also removing the interior cylinder from our model I was able to restore the butterfly modes as well. I can see why removing the constraints and forces would make a difference, but I don't understand how the interior region which was only used to aid in making the meshing in the central region more dense would have such a strong effect. With that removed as well, we get the following results using a finer mesh (plot attached):
ratioR Lowest Drumhead mode (Hz) Lowest Butterfly mode (Hz) Any lower modes (Hz)
0.1 8777.8 7253.3 3059.7
0.3 5636.9 7264.7 4202.3
0.5 6068.5 7308.0 5980.9
0.75 7341.5 7183.0 no
1 8110.2 5951.3 no
1.1 7302.8 5308.8 no
1.3 5744.3 4087.9 no
This data is more in line with what we would expect, and indicates that the optimal region is probably between ratioR=0.6 and ratioR=0.8. I will run the looping script we designed to map out how Umax changes when varying ratioR, particularly in this interesting region.
This line was used to generate the plot from data.
figure(1); plot(ratio,drum); hold on; plot(ratio,butterfly,'rs'); plot(ratio(1:3),others,'g'); hleg1=legend('Drumhead','Butterfly','Other'); title('Eigenfrequencies of mode types for different ratioR values'); xlabel('ratioR'); ylabel('lowest eigenfrequency'); shg;
Quote: 
Quote: 
I attempted to run the Matlab code built last week, but received an error that the mesh failed to generate on the inner domain/failed to respect edge element on geometry edge. The error occurred on the second computation with ratioR=0.2, but for ratioR=0.1 the simulation completed successfully. I reproduced the error in COMSOL and noted that for meshes with size at least as course as "course" or extremely fine, a mesh could be constructed, but for sizes in between an error was returned. This region of sizes which throw errors varies for different values of ratioR. By specifying an edge selection of the central line running down the frustum and another along the top edge of the inner cylinder and meshing points on those edges first the rest of the meshing is forced to use those nodes as a seed.
I also noted our COMSOL model had an error: ratioR was applied to the wrong face. This was corrected as was the .m file.
I added an eigenfrequency study and ran it for several values of ratioR, obtaining the following:
ratioR Lowest Drumhead mode (Hz) Lowest Butterfly mode (Hz)
0.1 3624.8 3802.4
0.3 4375.9 3785.4
0.5 4366.5 3650.0
.75 3709.0 3041.85
1 3000 2186.3
1.1 2736.2 1893.0
1.3 2247.6 1426.4
The mirror's front face is 0.17 [m]. Its Gaussian beam size is 0.0156 [m]. Its height is given by 0.01734[m^3]/(R^2*(1+ratioR+ratioR^2)) such that when the two faces have the same radii it is 0.2 [m] long.

When we talked today, I thought that these frequencies were a little low, so I ran a really simple COMSOL model of a cylinder with radius 0.17 m and height 0.2 m (attached).
I solved for the first 10 modes (the first 6 being trivial), and the frequencies I got are listed below:
Mode # Freq [Hz] Type
7 5950 Butterfly
8 5950 Butterfly
9 8106 Drumhead
10 8226 Shear
These are more in the frequency range that I expected.
**What material are you using? I used fused silica.
**I used 'finer' meshing. Perhaps your more complicated meshing is screwing with the solutions, although I'd expect it to make it stiffer, not softer.
**Are you applying force during the eigenvalue calculation? That might make it softer.
**Perhaps some of your constraints are moving the modes around. My model has zero constraints: it's just a cylinder floating in space.
Lets talk about this tomorrow.


Attachment 1: Eigenfrequency.png


52

Wed Jul 10 14:15:20 2013 
Emory Brown  Optics  General  Comparison of Umax and the frequency of the lowest eigenmode against ratioR 
I was able to run the looper program designed last week over the interval ratioR=[0.5, 0.9] with a step size of 0.01 and obtained the plot below with a minimal value of Umax=1.4723e10 J obtained for ratioR=0.71, which is about 3.2% less than the value of Umax=1.5202e10 J when ratioR=1. This is not very significant, but the shifting upward of the frequencies of the lowest eigenmodes as described in the previous post is. Additionally, these curves both have their optimal points in the same region, with a ratioR value of about 0.7.
I also ran a higher resolution test to see how the frequency of the lowest real eigenmode responded to changes in ratioR over this region. The plot generated is also attached with a maximal (most favorable) value of 7210.0 Hz at ratioR=0.74.
Data from both of these runs is included in the attached tarball stored as txt files delineated by tabs.
The next few things to try are substituting silicon for fused silica as the material in the model and seeing how the curves change and seeing the response from reversing the orientation of the optimal test mass (ratioR~0.72). I expect that the orientation change will result in much higher noise and a lower minimal real eigenfrequency. 
Attachment 1: UmaxValuesVsRatioR.png


Attachment 2: FrequencyLowestEigenmodes.png


Attachment 3: data.tar.gz

53

Wed Jul 10 17:52:18 2013 
Koji  General  General  How to force time steps for a time dependent simulation 
A common trap in COMSOL:
When you run a time series simulation, COMSOL tries to be "clever" by skipping some time steps
and returns you an interpolated result in order to accelerate the calculation.
But for certain types of applications, such as time series response analysis by applying periodic force/heat,
this "clever" adaptive algorithm messes up our calculation.
You can apply forced time steps by giving a specific option:
 In Model Builder, open your study. Then open your Solver
 Click "TimeDependent Solver" in Model Builder
 You get bunch of options in TimeDependent Solver pane.
Open "Time Stepping".
 Change "Method" from "Generalized alpha" to "BDF".
 Change "Steps taken by solver" to "Strict"
Then, you should be OK.
This can be done by
model.sol('sol1').feature('t1').set('tstepsbdf', 'strict');
in COMSOL with MATLAB. (Of course sol1 and t1 should be changed appropriately.) 
54

Thu Jul 11 19:16:12 2013 
Deep Chatterjee  General  General  TE noise problem  Matlab script 
The present aim is the calculation of the TE noise like Liu and Thorne. Although, the results of Liu and Thorne are for the case of an infinite mirror, in our case, we try
to model a test mass with large dimensions in an attempt to get closer to the result due to Liu and Thorne
Following my previous post  The presents results given by the COMSOL simulations gave the same profile as Liu and Thorne's result which has a frequency dependence
of 1/omega^2 but was displaced from the analytic result by a constant factor which was dependent on the applied pressure on the face of the test mass which should not be
the case.
Following the steps followed Liu and Thorne, we had constructed the test mass as a cylinder, large compared to the beam spot size. An oscillating pressure was applied
to one of the faces and the from the temperature gradient generated in the process due to strains, one can calculate the work dissipated. The process involves integrating
the gradient of temperature over the entire geometry and taking the time avrage. However, this did not give the correct results, so it was decide to extract the Fourier coefficient
of the signal and perform the integration on it, as done in the case of TR niose by Koji Arai.
I mention the steps it was done 
> The data extracted was stored in a 3D array in matlab. using mphinterp  each dimension for r, z, t
> The number of radial slices and z slices is defined by the user previously
> For each r,z value there was a time signal.
> The fourier coefficient corresponding to the pressure oscillation frequency was extracted from the last three cycles
(In theory only one cycle should suffice)
> The above step was done so that random data( which is small in magnitude) generated by COMSOL is avoided.
> The fourier coefficients are stored in a 2D array corresponding to r and z values.
> The integration was performed using trapz twice once on each dimension to get the total volume integral
> The rest of the calculation is same as the previous script  plugging in the prefactors in the formulae and plotting
**The primary problem with this script is that the data extraction takes significant amount of time  For a mirror with 200 times the radius of the beam spot and 200 radial slices, it takes close to 6 minutes
to evaluate work dissipated for each frequency. The script solves for 16 frequency values. The results for small number of radial slices does not follow the straight line profile. For larger number of slices
the program takes a longer time. The results have not been checked yet. 
Attachment 1: TE_model.m

% 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;
... 94 more lines ...

Attachment 2: almost_infinite_2D_axisym_TE_4.m

function out = almost_infinite_2D_axisym_TE_4(param)
% Putting all the parameters in the variables
r0 = param.beam.radius;
R = param.mirror.radius;
H = param.mirror.height;
dt = param.COMSOL.simulation_time_step;
fmod = param.COMSOL.fmod;
t_end = param.COMSOL.simulation_end_time;
p0 = param.constant.F0/(pi*r0^2);% The amplitude of Osc. pressure
T_amb = param.constant.T_amb;
... 121 more lines ...

55

Fri Jul 12 00:09:17 2013 
Emory Brown  Optics  General  Analysis of substrate Brownian noise in a silicon test mass with changing ratioR value 
I performed the same analysis we previously did on a test mass made of silicon (using the builtin silicon within COMSOL) and obtained the attached plots. The plots show more of a difference between the two solutions optimal points than for fused silica, but they are still quite close. All data is included in the attached tarball. The values for ratioR=1 or a cylindrical test mass place the frequency of the lowest real eigenmode at 8491 Hz and the Umax value of 6.5484*10^11 J. 
Attachment 1: SiliconEigenfrequencyPlot.png


Attachment 2: UmaxPlotSilicon.png


Attachment 3: Silicon.tar.gz

56

Fri Jul 12 00:31:11 2013 
Emory Brown  Optics  General  Reversed face mirror design 
We had noticed that when we ran an eigenfrequency analysis, there was often a node (see attached screenshot) on the face opposite the applied force. This indicated that the noise might be lowered by reversing the face the force was applied to, but it also seems like this would be equivalent to increasing the ratioR value and making the primary face smaller while retaining the mass of the test mass. It seems worthwhile to see how the Umax values compare for a beam incident on either face of a test mass in the optimal region with Umax=0.71. For the initial configuration, Umax=1.47228*10^10 J and changing the side of the mirror the force is applied to we get Umax=1.47148*10^10 J, which is about the same value. 
Attachment 1: node.png


57

Fri Jul 12 19:43:32 2013 
Deep Chatterjee  General  General  TE noise problem  Matlab script 
Quote: 
The present aim is the calculation of the TE noise like Liu and Thorne. Although, the results of Liu and Thorne are for the case of an infinite mirror, in our case, we try
to model a test mass with large dimensions in an attempt to get closer to the result due to Liu and Thorne
Following my previous post  The presents results given by the COMSOL simulations gave the same profile as Liu and Thorne's result which has a frequency dependence
of 1/omega^2 but was displaced from the analytic result by a constant factor which was dependent on the applied pressure on the face of the test mass which should not be
the case.
Following the steps followed Liu and Thorne, we had constructed the test mass as a cylinder, large compared to the beam spot size. An oscillating pressure was applied
to one of the faces and the from the temperature gradient generated in the process due to strains, one can calculate the work dissipated. The process involves integrating
the gradient of temperature over the entire geometry and taking the time avrage. However, this did not give the correct results, so it was decide to extract the Fourier coefficient
of the signal and perform the integration on it, as done in the case of TR niose by Koji Arai.
I mention the steps it was done 
> The data extracted was stored in a 3D array in matlab. using mphinterp  each dimension for r, z, t
> The number of radial slices and z slices is defined by the user previously
> For each r,z value there was a time signal.
> The fourier coefficient corresponding to the pressure oscillation frequency was extracted from the last three cycles
(In theory only one cycle should suffice)
> The above step was done so that random data( which is small in magnitude) generated by COMSOL is avoided.
> The fourier coefficients are stored in a 2D array corresponding to r and z values.
> The integration was performed using trapz twice once on each dimension to get the total volume integral
> The rest of the calculation is same as the previous script  plugging in the prefactors in the formulae and plotting
**The primary problem with this script is that the data extraction takes significant amount of time  For a mirror with 200 times the radius of the beam spot and 200 radial slices, it takes close to 6 minutes
to evaluate work dissipated for each frequency. The script solves for 16 frequency values. The results for small number of radial slices does not follow the straight line profile. For larger number of slices
the program takes a longer time. The results have not been checked yet.

I am replying to report of the modification made to the script from the last one. The problem of too many radial and z slices have been avoided by using a slicing based on a gaussian
distribution both along the r and z directions  i.e. the slicing is heavy in the region around the surface is applied and much thinner towards the edges of the cylinder. The slicing has
been separately handled by a function named giveSlices() which does the slicing and returns the values as two arrays corresponding to r and z.
The extent to which the fineness of slicing is controlled in the geometry is by controlling the parameter SD in the code which is the standard deviation of the gaussian according to which
the slicing is controlled as being fine throughout or fine in the centre and thin with increasing r.
This method reduces the number of slices by almost 2 orders for significant large values of the dimensions of the test mass while it is expected to give similar results since slicing the
test mass finely at the edge is not required as the gradient of temperature almost falls to zero.
It is suspected that the long time of simulation is attributed to calling mphinterp() command a large number of times. In my next modification, it will be tried to use the command just once
followed by the proper data extraction from the output given by the command.

Attachment 1: codes_TE_calc.zip

58

Sun Jul 14 15:13:46 2013 
Deep Chatterjee  General  General  TE noise problem  Matlab script 
Quote: 
Quote: 
The present aim is the calculation of the TE noise like Liu and Thorne. Although, the results of Liu and Thorne are for the case of an infinite mirror, in our case, we try
to model a test mass with large dimensions in an attempt to get closer to the result due to Liu and Thorne
Following my previous post  The presents results given by the COMSOL simulations gave the same profile as Liu and Thorne's result which has a frequency dependence
of 1/omega^2 but was displaced from the analytic result by a constant factor which was dependent on the applied pressure on the face of the test mass which should not be
the case.
Following the steps followed Liu and Thorne, we had constructed the test mass as a cylinder, large compared to the beam spot size. An oscillating pressure was applied
to one of the faces and the from the temperature gradient generated in the process due to strains, one can calculate the work dissipated. The process involves integrating
the gradient of temperature over the entire geometry and taking the time avrage. However, this did not give the correct results, so it was decide to extract the Fourier coefficient
of the signal and perform the integration on it, as done in the case of TR niose by Koji Arai.
I mention the steps it was done 
> The data extracted was stored in a 3D array in matlab. using mphinterp  each dimension for r, z, t
> The number of radial slices and z slices is defined by the user previously
> For each r,z value there was a time signal.
> The fourier coefficient corresponding to the pressure oscillation frequency was extracted from the last three cycles
(In theory only one cycle should suffice)
> The above step was done so that random data( which is small in magnitude) generated by COMSOL is avoided.
> The fourier coefficients are stored in a 2D array corresponding to r and z values.
> The integration was performed using trapz twice once on each dimension to get the total volume integral
> The rest of the calculation is same as the previous script  plugging in the prefactors in the formulae and plotting
**The primary problem with this script is that the data extraction takes significant amount of time  For a mirror with 200 times the radius of the beam spot and 200 radial slices, it takes close to 6 minutes
to evaluate work dissipated for each frequency. The script solves for 16 frequency values. The results for small number of radial slices does not follow the straight line profile. For larger number of slices
the program takes a longer time. The results have not been checked yet.

I am replying to report of the modification made to the script from the last one. The problem of too many radial and z slices have been avoided by using a slicing based on a gaussian
distribution both along the r and z directions  i.e. the slicing is heavy in the region around the surface is applied and much thinner towards the edges of the cylinder. The slicing has
been separately handled by a function named giveSlices() which does the slicing and returns the values as two arrays corresponding to r and z.
The extent to which the fineness of slicing is controlled in the geometry is by controlling the parameter SD in the code which is the standard deviation of the gaussian according to which
the slicing is controlled as being fine throughout or fine in the centre and thin with increasing r.
This method reduces the number of slices by almost 2 orders for significant large values of the dimensions of the test mass while it is expected to give similar results since slicing the
test mass finely at the edge is not required as the gradient of temperature almost falls to zero.
It is suspected that the long time of simulation is attributed to calling mphinterp() command a large number of times. In my next modification, it will be tried to use the command just once
followed by the proper data extraction from the output given by the command.

Following the previous reply regarding the use of mphinterp(), the appropriate changes in code were made. The overhead time due to the call to COMSOL was responsible for the longer
simulation time. Right now, for relatively large no. of total slices (about 27000) just the extraction of data happens within a few seconds. However, for the larger number of slices, the codes
will probably still take some time to complete. The codes are attached. 
Attachment 1: codes_TE_calc_(2).zip

59

Mon Jul 15 13:14:30 2013 
Emory Brown  Optics  General  Umax curves for varying height test masses 
I ran sweeps over the ratioR range from 0.5 to 0.9 for test masses with modified height, and thus mass, which were otherwise equivalent to our previous test mass. The attached plot was generated and indicates that the optimal ratioR value for the test mass was shifted significantly, for the half mass test mass into lower frequency regions, and for the 10x mass test mass into higher frequency ranges. I took the data collected from the half mass run and fit a polynomial of degree 2 to it in Matlab using the command p=polyfit(x,y,2) and obtained p=1.089*10^10*x^2  1.024*10^10*x +1.774*10^10 which has a minimum value for ratioR=0.47. This fits with the plot generated as it appears to be near its minumum value, but I should run the simulation more around this value to see if it fits the expected behaviour. I also did this fitting for the 10x mass curve, but based on the greater variance in the curve over this region of ratioR values I don't expect the predicted minimum to be as accurate, so if we want to find it we will have to search that region of ratioR values. For the 10x mass test mass, p=4.459*10^10*x^2  8.859*10^10*x + 6.289*10^10 which has a minimum of about 1.
I will run simulations in these regions overnight, and edit this post with the plots/results. 
Attachment 1: MassComparison.png


60

Mon Jul 15 16:05:29 2013 
Deep Chatterjee  General  General  TE noise problem  Matlab script 
Quote: 
Quote: 
Quote: 
The present aim is the calculation of the TE noise like Liu and Thorne. Although, the results of Liu and Thorne are for the case of an infinite mirror, in our case, we try
to model a test mass with large dimensions in an attempt to get closer to the result due to Liu and Thorne
Following my previous post  The presents results given by the COMSOL simulations gave the same profile as Liu and Thorne's result which has a frequency dependence
of 1/omega^2 but was displaced from the analytic result by a constant factor which was dependent on the applied pressure on the face of the test mass which should not be
the case.
Following the steps followed Liu and Thorne, we had constructed the test mass as a cylinder, large compared to the beam spot size. An oscillating pressure was applied
to one of the faces and the from the temperature gradient generated in the process due to strains, one can calculate the work dissipated. The process involves integrating
the gradient of temperature over the entire geometry and taking the time avrage. However, this did not give the correct results, so it was decide to extract the Fourier coefficient
of the signal and perform the integration on it, as done in the case of TR niose by Koji Arai.
I mention the steps it was done 
> The data extracted was stored in a 3D array in matlab. using mphinterp  each dimension for r, z, t
> The number of radial slices and z slices is defined by the user previously
> For each r,z value there was a time signal.
> The fourier coefficient corresponding to the pressure oscillation frequency was extracted from the last three cycles
(In theory only one cycle should suffice)
> The above step was done so that random data( which is small in magnitude) generated by COMSOL is avoided.
> The fourier coefficients are stored in a 2D array corresponding to r and z values.
> The integration was performed using trapz twice once on each dimension to get the total volume integral
> The rest of the calculation is same as the previous script  plugging in the prefactors in the formulae and plotting
**The primary problem with this script is that the data extraction takes significant amount of time  For a mirror with 200 times the radius of the beam spot and 200 radial slices, it takes close to 6 minutes
to evaluate work dissipated for each frequency. The script solves for 16 frequency values. The results for small number of radial slices does not follow the straight line profile. For larger number of slices
the program takes a longer time. The results have not been checked yet.

I am replying to report of the modification made to the script from the last one. The problem of too many radial and z slices have been avoided by using a slicing based on a gaussian
distribution both along the r and z directions  i.e. the slicing is heavy in the region around the surface is applied and much thinner towards the edges of the cylinder. The slicing has
been separately handled by a function named giveSlices() which does the slicing and returns the values as two arrays corresponding to r and z.
The extent to which the fineness of slicing is controlled in the geometry is by controlling the parameter SD in the code which is the standard deviation of the gaussian according to which
the slicing is controlled as being fine throughout or fine in the centre and thin with increasing r.
This method reduces the number of slices by almost 2 orders for significant large values of the dimensions of the test mass while it is expected to give similar results since slicing the
test mass finely at the edge is not required as the gradient of temperature almost falls to zero.
It is suspected that the long time of simulation is attributed to calling mphinterp() command a large number of times. In my next modification, it will be tried to use the command just once
followed by the proper data extraction from the output given by the command.

Following the previous reply regarding the use of mphinterp(), the appropriate changes in code were made. The overhead time due to the call to COMSOL was responsible for the longer
simulation time. Right now, for relatively large no. of total slices (about 27000) just the extraction of data happens within a few seconds. However, for the larger number of slices, the codes
will probably still take some time to complete. The codes are attached.

The dependence of the Power Spectral Density on the parameters F0, which is the pressure amplitude, has been somewhat corrected for in the codes attached. The mistake lay in the fact that
extraction of the Fourier coefficients was being done for the ht.gradTmag, the Fourier coefficients of which are not the correct ones to be used in this case since ht.gradTmag is not the correct time signal.
In the present case the Fourier coefficients are extracted from the r and z components, squared and added to carry out the calculation. However, the results are still off from the L&T results by order of 5.
The codes and the plots are attached.
The red values are the L&T result while the blue sketch is the simulation 
Attachment 1: codes_TE_calc.zip

61

Tue Jul 23 18:26:04 2013 
Deep Chatterjee  Optics  General  Comparison between Liu and Thorne Results and COMSOL results for TE noise 
In this post I report of the results of TE noise simulated by COMSOL for the TE noise of Infinite test masses.
The aim was to follow the procedure by Liu and Thorne in their analytic calculations so that the same model could be used for the other
geometries.
The simulation is done in a different way than the TR simulations. It was observed that the output given by COMSOL by the use of commands
like mphinterp() or taking an export resulted in certain discrepancies between the results computed in COMSOL and that read by MATLAB.
Thus, the volume integration of the temperature gradient is performed in COMSOL itself and the results of the integration for each time
are sent to files. Matlab read these values and time averages them to get the result as in the paper (Sec. 2 of Liu and Thorne).
The errors expected are
> Fourier analysis is not done at all. This would have involved exporting data which, as mentioned before is giving errors
> The numerical errors by COMSOL are therefore not filtered off.
> The plot differs from the analytic solution for larger frequencies over 3000 Hz.
> It is to be noted from the paper by Liu and Thorne that the TE noise for the finite and infinite case are not very different. In
fact the correction factor goes as O(1). Thus, differences between finite and inifinte cases are unlikely to be prominent
in the log scale plots
The codes are put as a zip file. Corrections made to the codes will be uploaded as a reply.

Attachment 2: TE_calc_Jul23.zip

Attachment 3: Jul_23.eps


67

Tue Jul 23 20:53:37 2013 
Deep Chatterjee  Optics  General  Comparison between Liu and Thorne Results and COMSOL results for TE noise 
Quote: 
In this post I report of the results of TE noise simulated by COMSOL for the TE noise of Infinite test masses.
The aim was to follow the procedure by Liu and Thorne in their analytic calculations so that the same model could be used for the other
geometries.
The simulation is done in a different way than the TR simulations. It was observed that the output given by COMSOL by the use of commands
like mphinterp() or taking an export resulted in certain discrepancies between the results computed in COMSOL and that read by MATLAB.
Thus, the volume integration of the temperature gradient is performed in COMSOL itself and the results of the integration for each time
are sent to files. Matlab read these values and time averages them to get the result as in the paper (Sec. 2 of Liu and Thorne).
The errors expected are
> Fourier analysis is not done at all. This would have involved exporting data which, as mentioned before is giving errors
> The numerical errors by COMSOL are therefore not filtered off.
> The plot differs from the analytic solution for larger frequencies over 3000 Hz.
> It is to be noted from the paper by Liu and Thorne that the TE noise for the finite and infinite case are not very different. In
fact the correction factor goes as O(1). Thus, differences between finite and inifinte cases are unlikely to be prominent
in the log scale plots
The codes are put as a zip file. Corrections made to the codes will be uploaded as a reply.

Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.
One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot.
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies. 
Attachment 1: Jul_23_bettermesh2.eps


66

Tue Jul 23 20:53:37 2013 
Deep Chatterjee  Optics  General  Comparison between Liu and Thorne Results and COMSOL results for TE noise 
Quote: 
In this post I report of the results of TE noise simulated by COMSOL for the TE noise of Infinite test masses.
The aim was to follow the procedure by Liu and Thorne in their analytic calculations so that the same model could be used for the other
geometries.
The simulation is done in a different way than the TR simulations. It was observed that the output given by COMSOL by the use of commands
like mphinterp() or taking an export resulted in certain discrepancies between the results computed in COMSOL and that read by MATLAB.
Thus, the volume integration of the temperature gradient is performed in COMSOL itself and the results of the integration for each time
are sent to files. Matlab read these values and time averages them to get the result as in the paper (Sec. 2 of Liu and Thorne).
The errors expected are
> Fourier analysis is not done at all. This would have involved exporting data which, as mentioned before is giving errors
> The numerical errors by COMSOL are therefore not filtered off.
> The plot differs from the analytic solution for larger frequencies over 3000 Hz.
> It is to be noted from the paper by Liu and Thorne that the TE noise for the finite and infinite case are not very different. In
fact the correction factor goes as O(1). Thus, differences between finite and inifinte cases are unlikely to be prominent
in the log scale plots
The codes are put as a zip file. Corrections made to the codes will be uploaded as a reply.

Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.
One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot.
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies. 
Attachment 1: Jul_23_bettermesh2.eps


65

Tue Jul 23 20:53:37 2013 
Deep Chatterjee  Optics  General  Comparison between Liu and Thorne Results and COMSOL results for TE noise 
Quote: 
In this post I report of the results of TE noise simulated by COMSOL for the TE noise of Infinite test masses.
The aim was to follow the procedure by Liu and Thorne in their analytic calculations so that the same model could be used for the other
geometries.
The simulation is done in a different way than the TR simulations. It was observed that the output given by COMSOL by the use of commands
like mphinterp() or taking an export resulted in certain discrepancies between the results computed in COMSOL and that read by MATLAB.
Thus, the volume integration of the temperature gradient is performed in COMSOL itself and the results of the integration for each time
are sent to files. Matlab read these values and time averages them to get the result as in the paper (Sec. 2 of Liu and Thorne).
The errors expected are
> Fourier analysis is not done at all. This would have involved exporting data which, as mentioned before is giving errors
> The numerical errors by COMSOL are therefore not filtered off.
> The plot differs from the analytic solution for larger frequencies over 3000 Hz.
> It is to be noted from the paper by Liu and Thorne that the TE noise for the finite and infinite case are not very different. In
fact the correction factor goes as O(1). Thus, differences between finite and inifinte cases are unlikely to be prominent
in the log scale plots
The codes are put as a zip file. Corrections made to the codes will be uploaded as a reply.

Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.
One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot.
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies. 
Attachment 1: Jul_23_bettermesh2.eps


64

Tue Jul 23 20:53:37 2013 
Deep Chatterjee  Optics  General  Comparison between Liu and Thorne Results and COMSOL results for TE noise 
Quote: 
In this post I report of the results of TE noise simulated by COMSOL for the TE noise of Infinite test masses.
The aim was to follow the procedure by Liu and Thorne in their analytic calculations so that the same model could be used for the other
geometries.
The simulation is done in a different way than the TR simulations. It was observed that the output given by COMSOL by the use of commands
like mphinterp() or taking an export resulted in certain discrepancies between the results computed in COMSOL and that read by MATLAB.
Thus, the volume integration of the temperature gradient is performed in COMSOL itself and the results of the integration for each time
are sent to files. Matlab read these values and time averages them to get the result as in the paper (Sec. 2 of Liu and Thorne).
The errors expected are
> Fourier analysis is not done at all. This would have involved exporting data which, as mentioned before is giving errors
> The numerical errors by COMSOL are therefore not filtered off.
> The plot differs from the analytic solution for larger frequencies over 3000 Hz.
> It is to be noted from the paper by Liu and Thorne that the TE noise for the finite and infinite case are not very different. In
fact the correction factor goes as O(1). Thus, differences between finite and inifinte cases are unlikely to be prominent
in the log scale plots
The codes are put as a zip file. Corrections made to the codes will be uploaded as a reply.

Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.
One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot.
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies. 
Attachment 1: Jul_23_bettermesh2.eps


63

Tue Jul 23 20:53:37 2013 
Deep Chatterjee  Optics  General  Comparison between Liu and Thorne Results and COMSOL results for TE noise 
Quote: 
In this post I report of the results of TE noise simulated by COMSOL for the TE noise of Infinite test masses.
The aim was to follow the procedure by Liu and Thorne in their analytic calculations so that the same model could be used for the other
geometries.
The simulation is done in a different way than the TR simulations. It was observed that the output given by COMSOL by the use of commands
like mphinterp() or taking an export resulted in certain discrepancies between the results computed in COMSOL and that read by MATLAB.
Thus, the volume integration of the temperature gradient is performed in COMSOL itself and the results of the integration for each time
are sent to files. Matlab read these values and time averages them to get the result as in the paper (Sec. 2 of Liu and Thorne).
The errors expected are
> Fourier analysis is not done at all. This would have involved exporting data which, as mentioned before is giving errors
> The numerical errors by COMSOL are therefore not filtered off.
> The plot differs from the analytic solution for larger frequencies over 3000 Hz.
> It is to be noted from the paper by Liu and Thorne that the TE noise for the finite and infinite case are not very different. In
fact the correction factor goes as O(1). Thus, differences between finite and inifinte cases are unlikely to be prominent
in the log scale plots
The codes are put as a zip file. Corrections made to the codes will be uploaded as a reply.

Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.
One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot.
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies. 
Attachment 1: Jul_23_bettermesh2.eps


62

Tue Jul 23 20:53:37 2013 
Deep Chatterjee  Optics  General  Comparison between Liu and Thorne Results and COMSOL results for TE noise 
Quote: 
In this post I report of the results of TE noise simulated by COMSOL for the TE noise of Infinite test masses.
The aim was to follow the procedure by Liu and Thorne in their analytic calculations so that the same model could be used for the other
geometries.
The simulation is done in a different way than the TR simulations. It was observed that the output given by COMSOL by the use of commands
like mphinterp() or taking an export resulted in certain discrepancies between the results computed in COMSOL and that read by MATLAB.
Thus, the volume integration of the temperature gradient is performed in COMSOL itself and the results of the integration for each time
are sent to files. Matlab read these values and time averages them to get the result as in the paper (Sec. 2 of Liu and Thorne).
The errors expected are
> Fourier analysis is not done at all. This would have involved exporting data which, as mentioned before is giving errors
> The numerical errors by COMSOL are therefore not filtered off.
> The plot differs from the analytic solution for larger frequencies over 3000 Hz.
> It is to be noted from the paper by Liu and Thorne that the TE noise for the finite and infinite case are not very different. In
fact the correction factor goes as O(1). Thus, differences between finite and inifinte cases are unlikely to be prominent
in the log scale plots
The codes are put as a zip file. Corrections made to the codes will be uploaded as a reply.

Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.
One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot.
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies. 
Attachment 1: Jul_23_bettermesh2.eps


68

Tue Jul 23 20:53:45 2013 
Deep Chatterjee  Optics  General  Comparison between Liu and Thorne Results and COMSOL results for TE noise 
Quote: 
In this post I report of the results of TE noise simulated by COMSOL for the TE noise of Infinite test masses.
The aim was to follow the procedure by Liu and Thorne in their analytic calculations so that the same model could be used for the other
geometries.
The simulation is done in a different way than the TR simulations. It was observed that the output given by COMSOL by the use of commands
like mphinterp() or taking an export resulted in certain discrepancies between the results computed in COMSOL and that read by MATLAB.
Thus, the volume integration of the temperature gradient is performed in COMSOL itself and the results of the integration for each time
are sent to files. Matlab read these values and time averages them to get the result as in the paper (Sec. 2 of Liu and Thorne).
The errors expected are
> Fourier analysis is not done at all. This would have involved exporting data which, as mentioned before is giving errors
> The numerical errors by COMSOL are therefore not filtered off.
> The plot differs from the analytic solution for larger frequencies over 3000 Hz.
> It is to be noted from the paper by Liu and Thorne that the TE noise for the finite and infinite case are not very different. In
fact the correction factor goes as O(1). Thus, differences between finite and inifinte cases are unlikely to be prominent
in the log scale plots
The codes are put as a zip file. Corrections made to the codes will be uploaded as a reply.

Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.
One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot.
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies. 
Attachment 1: Jul_23_bettermesh2.eps


69

Wed Jul 24 21:08:24 2013 
Deep Chatterjee  Optics  General  TR results for different dimensions 
In this post I simulate the procedure of calculating the TR noise for finite cavities as proposed by Heinert and check for a
match.
The technique of performing all necessary calculations in COMSOL and exporting the results was applied to the TR codes.
It was seen that the codes gives similar output as the technique of extraction of Fourier coefficients in place of time averaging
as has been done in the codes of Koji Arai. One can see the output as the present code runs to be similar to the previous ones
found in the SVN.
However, the results in the present case were off by a constant factor close to 100. This maybe due to some 'm'  'cm' or similar difference between
analytic calculations and COMSOL values of parameters. Although, it has not been found yet, the correction is hopeful to be
found soon.
The codes give results similar to the analytic result for other values of the mirror radius and beam radii (apart from the constant
factor I have mentioned above). One may have a look at the trend of the graphs between analytic and simulated values in the plots
attached. These plots are for the case when the mirror radius = 25m while the beam radius = 9 cm i.e. the original radii were 0.25m
and 9cm respectively i.e. the ratio has been changed by a order of 2.
As mentioned before the reason for the constant factor difference will be looked into. 
Attachment 1: Jul24.pdf


Attachment 2: Jul24_acalc.pdf


70

Thu Jul 25 11:36:17 2013 
rana  Optics  General  TR results for different dimensions 
PDF please instead of EPS or BMP or JFIF or TARGA or GIF or ascii art. 
71

Thu Jul 25 13:21:46 2013 
Deep Chatterjee  Optics  General  TR results for different dimensions 
Quote: 
In this post I simulate the procedure of calculating the TR noise for finite cavities as proposed by Heinert and check for a
match.
The technique of performing all necessary calculations in COMSOL and exporting the results was applied to the TR codes.
It was seen that the codes gives similar output as the technique of extraction of Fourier coefficients in place of time averaging
as has been done in the codes of Koji Arai. One can see the output as the present code runs to be similar to the previous ones
found in the SVN.
However, the results in the present case were off by a constant factor close to 100. This maybe due to some 'm'  'cm' or similar difference between
analytic calculations and COMSOL values of parameters. Although, it has not been found yet, the correction is hopeful to be
found soon.
The codes give results similar to the analytic result for other values of the mirror radius and beam radii (apart from the constant
factor I have mentioned above). One may have a look at the trend of the graphs between analytic and simulated values in the plots
attached. These plots are for the case when the mirror radius = 25m while the beam radius = 9 cm i.e. the original radii were 0.25m
and 9cm respectively i.e. the ratio has been changed by a order of 2.
As mentioned before the reason for the constant factor difference will be looked into.

The discrepancy related to the difference between the analytic and COMSOL results has been partially addressed. Attached is another
plot showing the comparison. The ratio this time between the COMSOL results and the analytic results is between 0.7  0.8. This difference
will be looked into. It is, however, observed that the difference is not a constant factor  it has to do with the model file. 
Attachment 1: Jul_25.pdf


72

Thu Jul 25 15:54:58 2013 
Deep Chatterjee  Optics  General  TR results for different dimensions 
Quote: 
Quote: 
In this post I simulate the procedure of calculating the TR noise for finite cavities as proposed by Heinert and check for a
match.
The technique of performing all necessary calculations in COMSOL and exporting the results was applied to the TR codes.
It was seen that the codes gives similar output as the technique of extraction of Fourier coefficients in place of time averaging
as has been done in the codes of Koji Arai. One can see the output as the present code runs to be similar to the previous ones
found in the SVN.
However, the results in the present case were off by a constant factor close to 100. This maybe due to some 'm'  'cm' or similar difference between
analytic calculations and COMSOL values of parameters. Although, it has not been found yet, the correction is hopeful to be
found soon.
The codes give results similar to the analytic result for other values of the mirror radius and beam radii (apart from the constant
factor I have mentioned above). One may have a look at the trend of the graphs between analytic and simulated values in the plots
attached. These plots are for the case when the mirror radius = 25m while the beam radius = 9 cm i.e. the original radii were 0.25m
and 9cm respectively i.e. the ratio has been changed by a order of 2.
As mentioned before the reason for the constant factor difference will be looked into.

The discrepancy related to the difference between the analytic and COMSOL results has been partially addressed. Attached is another
plot showing the comparison. The ratio this time between the COMSOL results and the analytic results is between 0.7  0.8. This difference
will be looked into. It is, however, observed that the difference is not a constant factor  it has to do with the model file.

The issue related to the difference between the analytic and simulated values has been resolved. The codes seems to give reasonable match
between the analytic and simulated case. There is, however, a difference between the formulas being used from the previous cases. Note that
the 1/2 in front of Eq.(15) of Heinert is a because the time average has already been considered. However, in the present codes, the volume
integral of grad_T is evaluated in COMSOL and exported as a function of time. It is then averaged in MATLAB. Thus the factor of 1/2 is to be
omitted in this case(see Liu and Thorne Eq.(5). The presence of this extra factor of 1/2 was giving error in the last upoaded plots. From the
relative difference plot, one can see the maximum difference between COMSOL and analytic results go upto 7% but for most of the graph
it is close to 1% which is a fair result. 
Attachment 1: Jul_25_1.pdf


Attachment 2: Jul_25_2.pdf


Attachment 3: relative_plot.pdf


73

Mon Jul 29 22:42:57 2013 
Deep Chatterjee  Optics  General  Avoiding transient solutions in the Computation of TE/TR noise 
An error was being encountered in the computation of the TR noise lately. It was observed that while running the simulations in the case of the materials which have a lower value of the thermal diffusivity (silicon / sapphire at room temperature), the simulated result were slightly off from the analytic result. On the other hand, if the simulation was run with a material of higher diffusivity(same materials at lower temperature), the match would be better. The reason being the transient solution not dying off significantly during the period of the simulation. Since a time average was being taken of the quantity integral{grad_T ^2}, the transient contributed to the integral. To get the correct value, the fourier coefficient of the time signal of integral{grad_T ^2}, was extracted at twice the frequency of the pressure oscillation. The reason being that the signal was squared. Extracting the response at this frequency after the integration is logical since the integration is over space while the response we extract is over time.
The same procedure was also applied to the TE noise calculation. However, this time we obtained similar result as the case where this procedure was not applied but a simple time averaging was performed. The tail of the plot, at high frequency, is still seen to deviate from the analytic result of Liu and Thorne as was the case previously. A plot is attached showing the spectrum for Fused silica at 290K. The conclusion being that the transients do not affect the TE noise calculation  the plot stayed the same even after filtering them out. This is probably because unlike the TR case which has a heat source present along the cylinder axis, the TE noise calculation involves applying pressure only on the face of the cylinder, and the transient do not contribute much to the volume integral.

Attachment 1: Jul29_2.pdf


74

Wed Jul 31 15:39:11 2013 
Deep Chatterjee  Optics  General  First try with paramter optimization for TE and TR noise profiles 
After the simulations have been found to match to a fair extent with the analytic results by Heinert and Liu and Thorne, the attempt is check out the parameters for which the TE and TR noise are
close to each other. This was done with the analytic results. The frequency range 10  1000 Hz was looked into. In between this range, the quantity that was minimized was the absolite value
of the logarithm of the ratio between the TR and TE noise. The fminsearch function was used to minimize the mentioned quantity. The parameters that were changed were  conductivity, thermo
optic coefficient and coefficient of linear expansion. The reason for choosing these three were 
> TR noise is independent of coefficient of linear expansion
> TE noise is independent of thermo optic coefficient
> The power law dependence on conductivity is different for the TE and TR cases as can be seen from the analytic expressions
once the code returned the optimized parameters, these values were plugged in and the results were plotted.
**Note that the minimization was done for frequencies between 10 to 1000 Hz 
Attachment 1: Jul31_param_opt.pdf


75

Thu Aug 8 17:17:19 2013 
Deep Chatterjee  Optics  General  Something like cancellation 
For the material parameters of Sapphire at 300K, the TE and TR Noise profiles, though not very close, lie close to within an order of magnitude. Sapphire has a positive coefficient of linear expansion. We just inverted the sign of this quantity
and the ran the codes that puts heat and pressure simultaneously to the test mass. The total noise looks to be lower than the TR noise which is greater.
Mirror radius 0.25[m]
Mirror height 0.46[m]
Beam Radius 0.09[m]
If we have physical parameters which make the two Noise sources come closer to each other and then flip the sign of alpha, we may be able to see some noise reduction to a greater extent. 
Attachment 1: suspected_cance_AUg8l.pdf


86

Tue Jun 24 14:35:42 2014 
Sam Moore  Optics  General  Trying to Verify the Heinert Model 

Attachment 1: 06_23_14.pdf


87

Tue Jun 24 17:05:24 2014 
Sam Moore  Optics  General  Trying to Verify the Heinert Model 
It does appear that the simplified model is only relevant for the simulations. To quote Heinert: "An efficient computation is only possible for the simple model, as the advanced model would require an element of size more than 10^{6 }." I have run Koji's code that replicates Heinert's figure 3. I have attached the resulting temperature distribution and noise amplitude curve. In the noise amplitude curve, the red line is the analytical result, while the dots are from COMSOL.
The next step is to convert this code to an efficient complex timeindependent solution. As stated before, my main concern here is whether COMSOL actually solves the right equation in the stationary case. 
Attachment 1: temperature.png


Attachment 2: noiseAmplitude_agreement.png


88

Sat Jun 28 21:59:11 2014 
Sam Moore  Optics  General  Difficulty with the COMSOL stationary module; Test Cases 
Here, I describe some test cases to see if COMSOL's solutions are agreeing with some simple analytical solutions. Right now, I have two plots showing COMSOL's solution and my analytical solution on separate plots. I will be plotting there difference to see if they really match up.

Attachment 1: 6_27_14.pdf


89

Sun Jun 29 15:37:18 2014 
Sam Moore  Optics  General  Difficulty with the COMSOL stationary module; Test Cases 
Quote: 
Here, I describe some test cases to see if COMSOL's solutions are agreeing with some simple analytical solutions. Right now, I have two plots showing COMSOL's solution and my analytical solution on separate plots. I will be plotting there difference to see if they really match up.

The following document shows the relative difference between these two plots. 
Attachment 1: 6_29_14.pdf


90

Sun Jun 29 20:25:44 2014 
Koji  Optics  General  Difficulty with the COMSOL stationary module; Test Cases 
What about this example? The result is easier to understand intuitively.
Consider a bar with the length of L.
Let's say there is no body heat applied, but the temperature of the bar at x=L is kept at T=0
and at x=0 is kept at T=T0 Exp[I w t].
The equation for the bar is
...(1)
Consider the solution with the form of T(x, t) = T(x) T0 Exp(I w t), where T(x) is the position dependent transfer function.
T(x) is a complex function.
Eq.1 is modified with T(x) as
With the boundary condition of
This can be analytically solved in the following form
where alpha is defined by
So kappa/Cp is the characteristic (angular) frequency of the system.
Here is the example plot for L=1 and alpha = 1 (red), 10 (yellow green), 100 (turquoise), 1000 (blue)
If the oscillation is slow enough, the temperature decay length is longer than the bar length and thus the temperature is linear to the position.
If the oscillation is fast, the decay length is significantly shorter than the bar length and the temperature dependence on the position is exponential.
Now what we need is to solve this in COMSOL 
93

Thu Jul 10 16:51:14 2014 
Sam Moore  Optics  General  Duan and Heinert Comparison 
(See Plots in attached document)
My plan has been to replicate Duan's numerical thermoconductive (TE + TR) phase noise plot presented in his paper (section V). I am trying to match Duan's analytical expression with Heinert's analytical expression. This requires some rescaling of Heinert's TR displacement noise. (I also needed to divide Heinert's expression by 4 pi to match the Fourier Transform convention. ) Duan's analytical expression for the phase noise is obtained by evaluating the triple integral given in equation 13 of the Duan paper "General Treatment of Thermal Noise in Optical Fibers".
It turns out that an additional factor of 2 multiplies the phase noise because Duan's Fourier Transform only takes into account positive frequencies; there are also negative frequencies that occur in equal amplitude.
This integral was evaluated in Mathematica due to numerical noise in MATLAB's calculation. The calculation in Mathematica was very slow, so the upper limits on the integral were truncated. The following plots in the attached document show the resulting noise profile agreements for two different upper limits.
If the residual for the highest upper limit is considered acceptable for a match between the two plots, then I will use Heinert's plot as a reference when using the COMSOL steadystate method for Duan's numerical case (Heinert's plot runs much faster). 
Attachment 1: 7_10_14.pdf


95

Mon Jul 14 19:09:14 2014 
Sam Moore  Optics  General  Using Heinert's Solution for Duan's Parameters 
I have plotted Heinert's analytical solution for TR noise using Duan's parameters. Since TO and TE noise can be found by simply rescaling TR noise, these have been included in the plot as well. The solid curve represents the analytical solution, while the tick marks represent COMSOL's solution. I have used COMSOL for both a 1D axisymmetric and a 3D model. Since Duan's cylinder has a radius of 125 microns, but a length of 1 m, the meshing was difficult for the 3D model. I ended up shortening the length of the cylinder, converting to the actual length when finally calculating the thermal noise. 
Attachment 1: oneD_duanParams_residualepsconvertedto.pdf


Attachment 2: threeD_duanParams_residualepsconvertedto.pdf


96

Mon Jul 14 19:14:31 2014 
Sam Moore  Optics  General  Duan and Heinert Comparison 
Quote: 
(See Plots in attached document)
My plan has been to replicate Duan's numerical thermoconductive (TE + TR) phase noise plot presented in his paper (section V). I am trying to match Duan's analytical expression with Heinert's analytical expression. This requires some rescaling of Heinert's TR displacement noise. (I also needed to divide Heinert's expression by 4 pi to match the Fourier Transform convention. ) Duan's analytical expression for the phase noise is obtained by evaluating the triple integral given in equation 13 of the Duan paper "General Treatment of Thermal Noise in Optical Fibers".
It turns out that an additional factor of 2 multiplies the phase noise because Duan's Fourier Transform only takes into account positive frequencies; there are also negative frequencies that occur in equal amplitude.
This integral was evaluated in Mathematica due to numerical noise in MATLAB's calculation. The calculation in Mathematica was very slow, so the upper limits on the integral were truncated. The following plots in the attached document show the resulting noise profile agreements for two different upper limits.
If the residual for the highest upper limit is considered acceptable for a match between the two plots, then I will use Heinert's plot as a reference when using the COMSOL steadystate method for Duan's numerical case (Heinert's plot runs much faster).

I have now plotted the DuanHeinert comparison for the case of an infinite upper bound. It turns out that the curves differ by a maximum of 10 percent for low frequencies. Such a discrepancy has been attributed to lack of experimental investigation into this regime (according to Duan). For our purposes, such a discrepancy is acceptable. We will therefore use the Heinert curve for subsequent calculations due to its faster computation time.

Attachment 1: duan_heinert_comparisonInfiniteepsconvertedto.pdf


97

Thu Jul 31 20:55:38 2014 
Sam Moore  Optics  General  Finding the Right Meshing for the TIR cavity 
In this document, I try to identify I good mesh by comparing the numerical solution from that mesh with my analytical model. Since there are problems with carrying out the analytical calculation, it is still not entirely clear which mesh should be used.

Attachment 1: 7_30_14.pdf


98

Sat Aug 2 00:22:34 2014 
Sam Moore  Optics  General  Finding the Right Meshing for the TIR cavity 
Quote: 
In this document, I try to identify I good mesh by comparing the numerical solution from that mesh with my analytical model. Since there are problems with carrying out the analytical calculation, it is still not entirely clear which mesh should be used.

I have refined the analytical calculation procedure, as outlined in this new document. The procedure indicates that the discrepancy between the analytical and numerical solutions are more likely attributed to meshing inaccuracies. 
Attachment 1: 8_1_14.pdf


108

Fri Jul 29 14:33:41 2016 
Joy Westland  Mechanics  General  A Tutorial in Importing SolidWorks Files and Applying a Gaussian Force in ANSYS for a LIGO Test Mass 
Here is a tutorial in importing SolidWorks into ANSYS and the steps needed to apply a Gaussian Force to the LIGO test mass that’s imported.
Using SolidWorks:
 Download the SolidWorks Zip Folders from the DCC
 https://dcc.ligo.org/login/index.shtml?entityID=https%3A%2F%2Fdcc.ligo.org%2Fshibbolethsp&return=https%3A%2F%2Fdcc.ligo.org%2FShibboleth.sso%2FLogin%3FSAMLDS%3D1%26target%3Dss%253Amem%253A15b6c314d87e3fa8b3768d89cb6b9836fe39c754 (LIGOD1000760v4)
 Open/export the zip folders into a different file
 For only the test mass and ears, open the file called: D0902456 ITM OPTIC WITH EARS ASSEMBLY. This will represent the bottom mass.
 An assembly can be imported into ANSYS but it’s easier to convert an assembly into a part before importing it into ANSYS. Also this will allow you to perform a split line.
 To do this, resave the assembly but before pressing save, change the file to PRT extension instead of Assembly.
 Open the PRT file of the test mass, make a new sketch on the face of the test mass.
 Draw a circle (Or whatever shape) on the test mass. In this case a 0.1 m circle was made in the middle of the test mass.
 From there go to “Insert” and then press “Curve” and then go to “Split Line”. Split line allows a user to project the sketch onto the part. This makes it possible to press the sketch separately from the part. This is a useful tool because in ANSYS, once split line has been used on a sketch, a user can press that shape separately on the face of the object. The importance of this tool allows a user to click on the circle separately and apply a force in that particular area.
 In “Split Line” press the “sketch” you want to project onto the “face” of the object.
Importing SolidWorks into ANSYS:
 Once the part has been made, choose the “static structural analysis”.
 In the geometry: Import the test mass that was converted into a part
 Once imported, press the geometry again and “edit”
 In the edit module, “generate” the part and once the part has loaded notice the circle can be pressed such that it is a different part of the test mass surface.
 Exit geometry module
Applying Gaussian Force:
 Enter the “model” module to edit the setup of the analysis
 Assign the material that you want for the test mass under “Geometry” à “Solid” under the “Material” and “Assignment”
 Under the “Coordinate Systems”, make a new coordinate system
 “Coordinate System:
 Click on the face, in this case the circle made on the test mass for the “geometry” and “apply”
 Under “Definition” à “Type” à Change from Cartesian to Cylindrical. This allows the Gaussian Force to be distributed correctly
 Under “Analysis Settings” Turn on “Large Deflection”
 Add “gravity” under “Static Structural”
 Treat the test mass is if it were the bottom mass so make sure gravity is pointing the correct way, such that the wires on the ears would extend upward.
 Insert “Fix” to the top surfaces of the ears
 Insert “Pressure”
 Press on the circle made to apply for the pressure
 Under “Definition” change the “Magnitude” to “Function” (the arrow at the end of the Magnitude entry bar)
 Once function has been activated: A Gaussian force of: 1/((3.141592654*0.0156^2)*2.718281828^((x/0.0156)^2)) was used in SI units. Note that ANSYS does not use symbols. Once that’s entered into the “Magnitude”, under “Function” the “Coordinate System” will appear. Change that to the coordinate system that was made in step 4 of the cylindrical system.
 Under “Graph Controls” Make sure that the XAxis is changed from Time to X.
 Pick a range for the graph and the number of segments that you want to look at.
 “Solve” the system and right click on “Solution” and evaluate the different results you need.

Attachment 1: Bottom_Test_Mass_With_Ears.PNG


Attachment 2: Applying_Parameters.PNG


Attachment 3: Applying_Gaussian_Force.PNG


114

Mon Jul 31 22:18:57 2017 
rana  General  General  using more than 12 cores in matlab 
Since 2014, the limit of 12 workers using the matlab parallel computing toolbox has been lifted. Today, I was able to get this to work. There's a trick.
Usually, when you start up matlab and run a parallel thing like 'parfor', it just uses a default profile 'local' which limits you to 12 workers. You can try to ask for more by doing 'parpool(40)' for 40 workers, but it will tell you that NumWorkers = 12 and you're out of luck. So instead:
myCluster = parcluster('local')
myCluster.NumWorkers = 40;
saveProfile(myCluster);
parpool('myCluster', 40)
It seems that it needs the max # of workers and the requested number of workers to be 40 to use 40, otherwise you'll just get 12 (as of matlab 2016a). 
Attachment 1: Screen_Shot_20170731_at_10.11.35_PM.png


127

Sat Mar 17 15:27:48 2018 
rana  General  General  file size >> small 
When saving your COMSOL files do these two things to make the files much smaller (good for saving in version control and sharing):
 File > Compact History
 Preferences > Files > Optimize for File Size (not speed)
